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Use  of  Long-Range  Weather  Forecasts  in  Ship  Routing 

G.  J.  HALTINER,  V.  E.  BLEICK  and  F.  D.  FAULKNER 

Naval  Postgraduate  School,  Monterey,  Calif. 

ABSTRACT 

Two  advances  in  the  calculus  of  variations  method  for  minimal- 
time  ship  routing  are  described.  The  first  is  a  scheme  for  con- 
structing ocean  wave  field  forecasts  which  may  be  expected  to 
have  considerable  skill  for  perhaps  8  days.  The  second  is  an  im- 
proved technique  for  varying  time-extremal  ship  tracks  toward 
admissibility.  Both  ideas  are  illustrated  by  calculating  the  op- 
timum track  ship  route  of  a  VC2AP3  vessel  on  a  trans-Pacific 
voyage.  Possible  future  developments  are  discussed. 
1 .  Introduction 

The  use  of  calculus  of  variations  methods  in  computing  minimal 
time  ship  routes  has  been  restricted  severely  in  the  past  by  the 
unavailability  of  ocean  wave  field  forecasts  for  extended  periods. 
In  an  initial  attempt  to  remedy  this  situation  it  is  proposed  here 
that  the  wave  forecasts  now  available  for  periods  up  to  two  days 
from  the  Fleet  Numerical  Weather  Facility  may  be  extrapolated  with 
considerable  skill  for  6  more  days  by  using  certain  5-day  and  30- 
day  forecasts  now  available  from  the  U.  S.  Weather  Bureau.  Bleick 
and  Faulkner  (1965)  gave  a  method  of  computing  minimal-time  ship 
routes  by  varying  time-extremal  ship  tracks  toward  the  admissi- 
bility of  reaching  a  desired  terminal  point.  Their  scheme  of  ex- 
tremal variation  is  refined  here  so  that  rapid  convergence  of  the 
track  iteration  process  is  assured.  An  example  of  wave  field  con- 
struction and  extremal  variation  is  given  for  a  Pacific  voyage. 

-1- 


2 .  Vave  Field  Construction 

The  scheme  of  incorporating  weather  forecasts  into  the  construc- 
tion of  a  10-member  computer-stored  time  series  of  wave  fields  for 
the  numerical  example  of  a  trans-Pacific  voyage  consists  of  the 
following  parts: 

a)  The  Fleet  Numerical  Weather  Facility  prepares  wave  analyses  at 
00Z  and  1 2Z  each  day,  as  well  as  operational  wave  predictions 
at  12-hour  intervals  for  periods  up  to  48  hours.  The  first  5 
members  of  the  time  series  consisted  of  the  analysis  at  1 2Z 

of  26  July  1966,  and  the  predictions  for  00Z  and  1 2Z  of  27  and 
28  July. 

b)  The  U.  S.  Veather  Bureau's  5-day  surface  pressure  forecast, 
issued  every  Monday,  Wednesday  and  Friday,  consisting  of  one 
sea-level  pressure  map  per  day  at  1230Z,  was  used  to  construct 
the  next  3  members  of  the  time  series.  The  last  3  maps  of  the 
forecast  issued  on  Wednesday,  27  July  1966,  were  used  to  de- 
termine surface  winds  and,  in  turn,  to  calculate  height,  per- 
iod and  direction  of  the  wind  waves  and  swell.  This  data  was 
used  for  time  series  members  at  1 2Z  of  29,  30  and  31  July. 

c)  The  U.  S,  Weather  Bureau's  30-day  forecast  was  utilized  for 
calculating  the  9th  and  10th  members  of  the  time  series.  Al- 
though not  published  for  use  outside  the  Weather  Bureau,  a 
copy  of  the  30-day  predicted  mean  sea-level  pressure  map, 
centered  at  the  middle  of  the  month  of  August  1966,  was  ob- 
tained. Surface  winds  were  estimated  from  this  single  map, 
and  again  wave  conditions  were  calculated.  The  calculations 
were  repeated  on  a  daily  basis  "using  the  same  map.  Since  the 

same  winds  are  used  repeatedly,  the  forecast  waves  reach  a 
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steady  state  within  a  few  days.  The  limited  amount  of  computer 
core  storage  made  it  necessary  to  terminate  the  time  series  by 
using  this  data  at  1 2Z  on  1  and  2  August.  The  last  member  of 
the  time  series  was  used  to  satisfy  any  later  need  for  wave 
fields . 
The  predicted  30-day  mean  pressure  chart  has  relatively  weak 
pressure  gradients,  as  would  be  expected  from  the  averaging  pro- 
cess. In  contrast,  the  individual  daily  charts  which  make  up  such 
a  mean  have  strong  gradients  in  general  ,  particularly  in  the  vi- 
cinity of  migratory  cyclones  or  low  pressure  areas.  These  systems 
have  strong  winds  and  high  seas  associated  with  them,  which  are 
reflected  in  the  30-day  mean  only  in  a  very  limited  fashion.  The 
forecast  procedure  did,  however,  show  considerable  skill  over 
the  use  of  long  term  monthly  mean  charts.  Nevertheless  it  is  de- 
sireable  to  seek  additional  ways,  possibly  more  accurate,  of  pro- 
viding wave  estimates  for  the  latter  part  of  a  voyage  extending 
beyond  5  days.  Possible  future  developments  of  this  kind  are 
discussed  later. 
3 .  Variation  of  Extremals 

As  in  the  previous  work  of  Bleick  and  Faulkner  (1965)  let  the 
differential  equations  of  a  ships  motion  in  the  stereographic 

plane  be 

x=V  cosp,   y=V  sinp.  (1) 

It  was  shown  that  the  ship  track  direction  angle  p  on  a  time- 
extremal  route  is 

p=arctan(|j/x) +arctan(V  /V)  .         (2) 

Here 

X=X-|  cosa+XpSina  (3) 

and 

1-1=1-1-,  cosa+ii-sina  (4) 
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arc  linear  combinations  of  the  linearly  independent  solutions 
X.,u1  and  Xp,Mp  of  the  adjoint  system 

X+V  (Xcosp+usinp)=0  (5) 

U+V  (Xcosp+|jsinp)=0 ,  (6) 

approximately 
and  ct  is, the  departure  angle  between  the  ship  track  and  the  Ox 

coordinate  axis  at  the  t=0  initial  point  of  the  voyage.  In  the 
previous  work  on  varying  a  time-extremal  ship  track  toward  the 
admissibility  of  reaching  a  desired  terminal  point  the  variation 
6p  was  considered  to  be  dependent  on  the  variation  6a  only.  This 
is  a  convenient  approximation  to  avoid  mathematical  complica- 
tions, but  its  use  may  lead  to  a  marked  slowing  down  of  the 
Newton-Raphson  track  iteration  process.  If  this  approximation  is 
abandoned  in  computing  the  variations  6x,6y  of  a  time-extremal 
ship  track  solution  of  (1)  and  (2),  then  the  dependence  of  6p  on 
all  of  the  variations  6x,  6y ,  6X-,  ,  ^X^,  6fi1  ,  5\i~    and  &a  must  be 
considered.  This  complete  variation  6p  is  found  from  (2)  to  be 

6p=[(V  /V)  6x+(V  /V)  6y+S2(x6C-u6e+E6a)]/D         (7) 

where 

6^=6X1 cosa+6x2sina,  (8) 

5c=6u1 cosa+6u2sina,  (9) 

s2=[i+(vp/v)2]/U2+M2),  (10) 

D  =  1+(V  /V)2-(V  /V)  (11) 

P        P    P 

E  =  X1M2-X2u1 .  (12) 

The  variation  of  (1),  the  time  differentiation  of  (8)  and  (9), 
and  use  of  (5),  (6)  and  (7)  give  the  following  non-homogeneous 
system  of  equations  to  solve  for  the  desired  variations  6x,6y: 
6x=[V  cosp-uQW  ]6x+[V  cosp-nQV  ]6y+uS2[Q(u6?-x6c) -E6a] ,   (13) 

A      x     y      y. 

6y=[v  sinp+XQV  ]6x+[V  sinp+XQV  ]6y-XS2[Q(u6?-X6c)-E6cc],   (14) 

A      x     y      y 
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-6£=S-1[V   +VD_,V  2]6x  +  S"1[V   +VD_,V  V  ]6y 
xx       x  xy       x  y  J 

+  [V  cosp-ugv  ]6p  +  [V  sinp+\QW  16C  +  QEW  6a ,         (15) 

.X.  A.  Jv  J\.  J\. 

-6C=S_1[V      +VD_1V  W    ]6x    +   S~1[V      +VD~1W   2]6y 
yx  y  xJ  yy  y   J  J 

+   [V   cosp-uQW    ]6?    +   [V    sinp+XQWr]6C   +  QEW   6a,  (16) 

*y  \j  *y  %j  \/ 

where  Q=VS/D  and  V=V  /V.  Equations  (13)  to  (16)  are  integrated, 
with  zero  initial  values  at  the  t=0  initial  point  of  the  track 
and  with  6a=1 ,  to  obtain  the  variations  6x(T)  and  6y(T)  at  the 
t=T  terminal  point.  These  variations  are  really  the  partial  de- 
rivatives dx(T)/da  and  dy(T)/da  since  we  have  taken  6a=1  .  The 
Newton-Raphson  equations  for  determining  AT  and  6a  on  a  varied 
time-extremal  track,  which  attempt  to  reduce  the  terminal  errors 
<-ix(T)  and  uy(T)  of  the  previous  extremal  track,  are  then 

x(T)aT+[Sx(T)/9a]6a=Ax(T) 

(17) 
y(T)AT+[3y(T)/9a]6a=ay(T) . 

In  the  numerical  integration  of  (1),  (5),  (6)  and  (13)  to  (16) 

use 
it  is  desir  able  to,  a  wave  field  interpolation  formula  which  will 

guarantee  as  far  as  possible  the  continuity  of  aJ 1  terms  of  these 
equations  where  any  of  x,y,t  assume  grid  values.  A  method  which 
achieves  this  when  interpolating  in  the  time  dimension  was  given 
by  Bleick  and  Faulkner  (1965).  The  method  given  there  for  inter- 
polating in  the  grid  of  the  Oxy  stereographic  plane  will  not  give 

the  desired  continuity  oi  V   ,  V    and  V    of  (15)  and  (16).  The 

J      xx    xy       yy 

16-point  interpolation  formula  used  here  to  guarantee  the  con- 
tinuity of  V  and  all  its  first  and  second  order  partial  deriva- 
tives with  respect  to  x  and  y,  except  V   ,  is  obtained  from  the 

xy 

4>'4  matrix  F(x,y),  whose  four  rows  and  columns  of  function  entries 
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correspond  to  four  successive  x  and  y  grid  values  respectively. 
The  interpolation  mesh  cell  is  the  central  cell  of  the  array, 
with  x  and  y  measured  from  the  cell  center,  and  with  the  mesh 
distance  considered  to  be  two  units.  The  formula  is 

F(x,y)  =  P(x)  F  P'(y)/1024  (18) 

where  the  row  matrix  P(x)=[P1  ,P?  ,P.  ,P . ]  has  the  elements 

P1=(x2-1)(x-1)2(x+2), 

P   =( 1 -x) (3x4+3x3-9x2-7x+1 8) , 
2  (19) 

P3=(x+1 ) (3x4-3x3-9x2+7x+18) , 

P4=(x2-1)(x+1)2(2-x), 

and  the  prime  indicates  matrix  transposition.  This  matrix  type  of 
interpolation  gave  excellent  results  in  the  numerical  example 
which  follows  despite  its  inability  to  give  continuity  of  V    at 

xy 

grid  values  of  x  or  y.  In  contrast  with  the  earlier  work  it  was 
found  desir  able  to  evaluate  the  various  derivatives  of  V=mv  by 
explicit  differentiation  of  the  solution  of  the  quadratic  equa- 
tion in  v  for  the  elliptical  polar  velocity  diagram. 
4 .  Numerical  Example 

Figure  1  illustrates  the  result  of  the  new  methods  of  wave 
field  construction  and  time-extremal  track  variation  in  the  case 
of  a  trans-Pacific  voyage  of  a  VC2AP3  vessel.  The  elliptical  pol- 
ar velocity  diagram  used  was  based  upon  the  work  of  James  (1959). 
The  minimal-time  track  starts  from  154E,  41N  at  1 200Z  on  26  July 
1966  and  ends  at  123V,  38N  at  0828Z  on  4  August,  with  circles  in- 
dicating successive  positions  of  the  vessel  at  8-hour  intervals. 
The  solid  line  is  the  great  circle  route  obtained  by  integrating 
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( I )  using 

cosp=n[(31-J)/31  .205]+m 

(20) 

sinp=n[(I-3l)/31 ,205]-l 
where  I, J  are  the  Fleet  Numerical  Weather  Facility  stereographic 
plane  grid  indices,  and  l,m,n  are  the  direction  cosines  of  the 
normal  vector  to  the  great  circle  plane.  These  cosines  are  com- 
puted from  the  normalized  cross  product  of  the  vectors  from  the 
earth's  center  to  the  initial  and  terminal  points  of  the  track. 
Because  of  the  rather  calm  prevailing  seas  there  was  no  signifi- 
cant time  difference  between  the  geodesic  and  minimal-time  routes, 
but  the  latter  did  show  a  reduction  in  the  wave  heights  encount- 
ered as  indicated  in  Fig.  1.  The  example  illustrates  the  advant- 
age of  the  new  method  of  extremal  variation  in  that  the  Newton- 
Raphson  Eqs.  (17)  were  used  as  they  stand  without  convergence 
difficulties,  i.e.  without  resorting  to  the  delayed  approach  to 
the  limit  scheme  of  using  only  some  fraction  of  6a  on  the  next 
track  iteration.  The  new  method  also  permitted  the  use  of  a  rather 
large  4-hour  time  step  in  the  numerical  integrations,  with  con- 
sequent gain  in  the  speed  of  the  track  iteration  process.  The 
Fortran  computer  program  may  be  obtained  from  the  authors. 
5 .    Concluding  Remarks 

Another  possibility  for  predicting  wave  fields  for  extended  per- 
iods, which  appears  to  have  promise,  is  to  utilize  a  wave  clima- 
tology. This  could  consist  of  utilizing  the  wave  analyses  now  be- 
ing prepared  daily  at  the  Fleet  Numerical  Veather  Facility  to  com- 
pute mean  wave  height,  direction  and  period  as  a  function  of  lati- 
tude and  longitude  for  each  month  of  the  year.  These  data  could 
then  be  compared  with  those  derived  from  the  Veather  Bureau  30- 
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day  sea  Level  pressure  forecasts  in  order  to  ascertain  the  best 
source  of  wave  data  for  trans-oceanic  ship  routing.  Such  a  wave 
climatology  would  have  other  applications  in  shipping  operations. 
A  further  refinement  in  the  development  of  such  a  wave  climatol- 
ogy might  consist  of  the  preparation  of  mean  wave  characteristics 
not  only  as  a  function  of  latitude,  longitude  and  month,  but  also 
separated  according  to  weather  type.  The  latter  are  determined 
largely  according  to  the  main  storm  tracks  which  vary  from  week 
to  week  as  well  as  with  season.  Such  a  climatology  would  obvious- 
ly take  more  effort  to  prepare,  but  would  be  a  very  valuable  aid 
in  ship  routing. 

Finally,  it  should  be  mentioned  that  a  number  of  groups  are  ex- 
perimenting with  long-range  weather  prediction  by  numerical  inte- 
gration of  the  hydrodynamical  equations.  It  is  expected  that 
eventually  such  predictions  will  show  skill  for  perhaps  several 
weeks,  and  thus  day-by-day  wave  forecasts  could  be  made  available 
for  the  entire  period  of  a  trans-oceanic  voyage. 
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8.  Appendix: 

Computer  Program  Input  and  Output 
The  program  was  written  for  the  CDC  1604  computer  in  Fortran 
I  963.  The  magnetic  tape  input  to  the  main  program  VC2AP3  of  this 
appendix  is  from  logical  unit  1.  This  input  tape  was  prepared  by 
program  TAPE  of  this  appendix,  and  involves  conversion  of  fixed- 
point  data  to  floating-point  data  by  the  normalizing  operation 
(addition  to  floating-point  zero)  of  the  CDC  1604.  The  first 
floating-point  BCD  record  of  the  tape  input  to  VC2AP3  contains 
the  data  required  to  plot  the  stereographic  plane  map  grid  of 
lines  of  longitude  and  circles  of  latitude  by  calling  subroutine 
DRAW ,  described  in  Naval  Postgraduate  School  Technical  Report/ 
Research  Paper  No.  73.  The  second  and  last  floating-point  BCD 
record  of  the  tape  input  to  VC2AP3  contains  the  three  18x32x10 
wave  field  time-series  arrays  XHT ,  CSK  and  SNK  corresponding  to 
the  wave  height  H,  and  the  wave  direction  cosines  cosK  and  sines 
sinK  described  in  the  first  reference.  The  dimensions  18  and  32 
correspond  to  the  FN¥F  stereographic  plane  grid  point  indices  of 

8<I<25  in  the  direction  of  the  1 0E  meridian  and  1 6<J<47  in  the 

10 
direction  of  the  1 00E  meridian.  The  dimension .  corresponds  to  the 

time  series  of  wave  fields  described  in  Section  2.  The  VC2AP3 
program  will  not  work  unless  all  points  of  a  ship  route,  includ- 
the  initial  and  terminal  points,  are  within  a  smaller  16x30  rect- 
angle defined  by  9<K24  and  17<J<46.  A  local  coordinate  system  is 
set  up  with  the  origin  0  at  1=7  and  J=1 5 ,  with  the  Ox  and  Oy  axes 
in  the  direction  of  increasing  I  and  J  respectively.  The  smallest 
values  of  x  and  y,  corresponding  to  1=8  and  J=1 6 ,  are  therefore 
x=1  and  y=1 .  The  punched  card  input  to  VC2AP3 ,  immediately  after 
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the  EXECUTER  control  card,  contains  the  following  data: 

Card  I:  First  line  of  the  TI=IT  title  in  format  (6A8/)  for  the 
map  produced  by  subroutine  DRAW  in  Statement  1 1 .  See 
example1  in  the  Fortran  card  listing  of  this  appendix. 

Card2:   Format  (8A8fI3)  of  which  6A8  is  the  second  line  of  the 

map  title  TI .  The  remaining  part  of  the  format  is  A8  for 
the  DATE=KATE  of  the  routing  computation,  A8  for  eight 
blank  Hollerith  characters  for  a  null  label  AL=LA  on  the 
map  grid  plot,  and  13  for  the  NST  total  number  of  ships 
to  be  routed.  The  DATE  of  the  routing  computation  cor- 
responds to  the  1 2Z  hour  of  the  first  member  of  the  time 
series  described  in  Section  2.  See  example  in  the  For- 
tran card  listing  of  this  appendix. 

Following  these  two  input  cards  there  a  groups  of  either  6  cards 
or  one  card  for  each  ship  routed  by  VC2AP3 ,  depending  on  whether 
or  not  the  option  to  plot  an  earlier  route  of  a  particular  ship 
is  elected . 

Card  3:  Format  (A4 , 2A8 ,F3 . 0 ,F6 . 1  ,F5 . 1  ,F6 . 1  ,F5 . 1  ,F3 .0 , 11  , 212 ,11  , 
F8.5,F6.3)  with  example  in  the  Fortran  card  listing  of 
this  appendix.  The  first  A4  is  the  GL=LG  ship  identifi- 
cation number  with  column  1  of  the  card  blank,  used  by 
subroutine  DRAV  to  label  the  terminal  point  of  a  ship 
route.  The  first  A8  is  the  DATEX=KATEX  date  on  which  the 
ship  leaves  the  initial  point  of  its  route.  The  second 
A8  is  the  FL=LF  Julian  date  of  departure,  with  blanks  in 
columns  17  to  20  inclusive,  used  by  subroutine  DRAV  to 
label  the  initial  point  of  a  ship  route.  The  F3.0  is  the 
HR  hour  of  ship  departure  from  its  initial  point  measur- 
ed from  1 2Z  on  the  DATE  of  routing,  i.e.  from  the  1 2Z 
hour  of  hour  of  the  first  member  of  the  wave  field  time 
series.  The  F6.1,  F5.1,  F6.1,  F5.1  formats  are  the  long- 
itudes and  latitudes  of  the  initial  and  terminal  points 
of  theroute,  XLG1 ,  XLT1 ,  XLG2 ,  XLT2,  with  longitudes 
considered  positive  if  east  of  the  Greenwich  meridian. 
The  F3 . 0  format  provides  for  the  RMUL  convergence  factor 
which  can  be  used,  to  divide  6a  by  RMUL  before  accepting 
it  for  the  next  track  iteration.  The  Fortran  card  list 
of  this  appendix  shows  RMUL=1  indicating  that  there  are 
no  convergence  difficulties  in  the  present  revision  of 
VC2AP3.  The  11  format  provides  for  NN  which  is  either  1 
or  zero  according  as  the  option  to  plot  an  earlier  route 
of  the  ship  is  or  is  not  elected.  The  first  of  the  212 
format  is  for  the  NSTEP  reciprocal  of  the  time  step  used 
in  the  integration  process,  with  6  of  the  card  listing 
indicating  a  step  of  1 /6  of  a  day.  The  second  12  format 
specifies  the  number  LMAX  of  iterations  allowed  in  det- 
ermining the  ship  route.  The  remaining  formats  11,F8.5, 
F6.3  provide  for  NP ,  PALF  and  PT  described  in  NPS  Tech. 
Report/Res .Paper  73,  but  found  unnecessary  in  present 
revision  of  VC2AP3  with  consequent  blanks  in  the  card 
listing. 

Card  3  is  followed  by  5  cards  punched  out  by  the  statements  on 
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cards  3  37  to  341  of  an  earlier  use  of  VC2AP3  if  the  option  NN=1 
to  plot  an  earlier  track  of  the  ship  has  been  elected.  If  NN=0 
on  card  3 ,  the  remaining  data  cards  of  the  input  deck  refer  to 
other  ships  to  be  routed. 

The  output  of  VC2AP3  contains  a  map  grid  for  each  vessel  rout- 
ed, shown  in  Fig.1,  produced  by  the  CALL  DRAV  of  statement  11.  If 
the  option  NN=1  has  been  elected,  statements  16  to  18  plot  an  ear- 
lier route  of  the  ship  using  plus  signs  for  daily  positions  and 
an  identifying  Julian  day  mark  for  the  initial  point.  Statements 

19  to  44  cause  a  geodesic  route  to  be  computed  and  plotted  as  a 

approximate 
solid  line  as  shown  in  Fig.1  where  the  initial  angle  of  departure 

ALF=a ,  measured  from  the  Ox  axis,  is  indicated  also.  The  terminal 
point  of  the  geodesic  route  is  marked  by  the  GL=LG  identification 
number  of  the  ship.  One  purpose  of  the  geodesic  route  computation 
is  to  find  first  approximations  to  the  time  T  and  departure  angle 
ALF  used  in  the  LMAX  iterations  toward  a  minimal-time  track  of 
statements  45  to  69.  Another  purpose  is  to  provide  a  standard  of 
comparison  for  the  effectiveness  of  the  minimal-time  routing.  The 
geodesic  route  computation  is  abandoned  if  any  point  of  the  route 
falls  outside  of  the  rectangle  9<K24  and  1 7<J<46 ,  but  the  route 
within  this  rectangle  is  plotted  on  the  map  grid.  The  minimal- 
time  route  computation  is  initiated  by  statement  45  only  if  the 
entire  geodesic  route  has  been  computed  successfully.  The  format 
of  statement  71  is  printed  if  the  LMAX  iterations  result  in  a 
terminal  point  more  than  100  nautical  miles  from  the  desired  de- 
stination, together  with  advice  about  how  to  improve  convergence. 
Experience  to  date  on  trans-Pacific  routes  indicates  that  it  is 
desireable  to  use  LMAX=1 0  and  RMUL=1  .  The  tabulated  daily  posi- 
tion, wave  height  and  direction  for  the  last  of  the  LMAX  itera- 
tions are  printed  under  the  format  of  statements  73  and  74,  with 
example  on  the  following  page.  Five  cards,  shown  on  the  following 
page,  are  punched  under  the  formats  of  statements  76  and  78, 
which  may  be  used  for  some  later  plot  of  the  track  if  the  NN=1 
option  of  card  41  of  VC2AP3  is  elected  in  a  later  routing  of  the 
ship.  The  statement  of  card  342  of  VC2AP3  causes  the  daily  track 
positions  to  be  plotted  on  the  map  grid  as  in  Fig.1 ,  with  the  LF 
=FL  Julian  day  identification  of  the  initial  point.  Statement  80 
continues  the  M=1 ,NST  loop  for  the  routing  of  the  next  ship. 
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PRINT  output  for  J2Q7  route  of  ship  666 

TCTAL  NU^eER  CF  VC2AP3  SHIPS  ROUTED  =    1 

CN  JUL26,66  FROP  12Z  TC  24Zf  *ND  CN  FOLLOWING  OAY  FRCP  OOZ  TC  III 

PCUTE  CF  SHIP   666  BEGINS  CN  JUL26,66,  JULIAN  CATE  =  J207, 
C  HCLRS  /JFTER  12C0Z  CN  JUL26,66 
FRCf  LONGITUDE  =   154. G  ANC  LATITUCE  =    41. C 
TC  LONGITUDE  =  -123.0  AND  LATITUCE  =    38. C 

R*LL  =     1    LPAX=1C    NSTEP=  6    NN=0    NP  =  C 


L     M       AL  F        T        X(N1)     XFIN      V(M)     YFIN 

C     55    -1.72*71     8.856    13.62C    13.620     A. 669     4.869 

1  55    -1.72471     8.856    14.112    13.62C     5.C20     4.869 

-1 .78  ICC     8.855 

2  55    -1  .78  ICC     8.855    13.915    13.620     4.96C     4.869 

-2.30619     8.85C 

3  49    -2.30619     8. COG     2.343    13.62C     7.467     4.869 

-1.49355     7.587 

4  47    -1.49355     7.587    16.95C    13.62C     9.125     4.869 

-2.29251     8.649 

5  49    -2.29251     8.CCC     2.516    13.62C     7.4C8     4.869 

-1.4263  9     7.614 

6  45    -1.42639     7.167    16.915    13.620    1C.C77     4.869 

—  2  C  1 C  Q  c     8  787 

7  52    -2.'01CC5     8!50C     8.146    13.62C     5.CC7     4.869 

-1.8C262     8.497 

8  52    -1.8C262     8.497    13.544    13.620     5.713     4.869 

-1.P2C49     8.852 

9  55    -1.82C49     8.852    13.5C1    13.62C     4.831     4.669 

—  1  P  1  *■  P  ^-      fl  P  *5  "^ 

10  55         -U81506  81853  13.623  13.620  4.669  4.869 

-1.81522  8.853 

CAYS        LCNCI     LATI-      WAVE      WAVE  CIRECTICN 
CF  TRAVEL     -TUDE     TLDE      HEIGHT       FRCP  NCRTH 

0  154. C  41. C  6.3  9  221 

l.CO  162.9  43.3  4.93  224 

2. CO  172.4  45. C  4.43  166 

3. CO  -177.6  46. C  6.64  2C2 

4. CO  -167.6  46.2  9.52  257 

5. CO  -157.6  45.8  9.81  251 

6.00  -148.  C  44.9  4.92  234 

7. CO  -138.6  43.3  3.93  168 

8. CO  -129.9  40.7  5.14  336 

8.e5  -123. C  38. C  2.81  331 

GPAPF  TITLEC 

JCe  C574         ELEICK  BO  6 

VC2AP3  APRIL  14         1967 

FAS  BEEN  PLCTTEC. 

PUNCH  output  for  J207  day  route  of  ship  666 

666     J207     10 

10.1210.1110.2510.5410.9711.5812.3913.3414.3415.30  666  J207  1 

0000000000  666  J207  2 

21.9719.6817.4415.2513.1010.99  8.93  6.91  4.88  3.15J207  666  1 

000000000            0J207  666  2 
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-C00P.B0X6.  BLEICK* I/l/0/49/S/56/57/E/45=54*15»10000*0,  VC2AP3  -  14  APR  67. 

-BINARY, 56. 

(RELOCOM. 

-FTN*L*A,E. 

PROGRAM  VC2AP3 

C   YVARS( 1 )=LAMBDA1   YVARS(2)=MU1     YVARS ( 3 ) =LAMBDA2     YVARS(4)=MU2 

C   YVARS(5)=X         YVARS(6)=Y       YVARS(7)=S  OR  VARX   YVARS ( 8 )= VARY 

C   YVARS(9)=VARXI    YVARS (  10  )  =VARZET A 

DIMENSION  X(900) ,Y<900)*RX(10,90),RY(10,90),IT(12).TI(12)»C(4),  1 

+  AK(4,10>  »DY( 10)  ,D(20)  2 

COMMON     YC( 10) *LR» A  * B *CC ,H »CK *SK  3 

COMMON/L1/XHK5760)  ,  CSM5  760)  ,  SNK(  5760)  ,  TC  4 

COMMON/ L5 /COST »S I  NT , V ,CAPV ,CAPVX *CAPVY , CAPVXX ♦ CAPVX Y *CAPVYY »  5 

+  VPBVX,VPBVY,DIV*RBV,CAPVP  6 

COMMON/L6/YVARS< 10) »XK*XLG*XLT  7 

EQUIVALENCE  (  I  T  ,T  I  ) , ( LA »AL  )  » (KATE* DATE) , ( LP , PL ) > ( LG *GL ) ♦ ( LF » FL )  ♦  8 

+  (X,RX ) ♦( Y*RY) ♦ (KATEX»DATEX)  9 

REWIND  1  10 

C( 1 )  =  0.0  11 

C(2)  =  0.5  12 

C(3)  =  0.5  13 

C(4)  =  1.0  14 

C  READ  MAP  GRID  DATA  FOR  DRAW  SUBROUTINE*  AND  WAVE  FIELD  ARRAYS 

REAO(lt3)  <X(  I  ) ,1  =  1*390) . (Y< I )  ,1  =  1*390)  15 

3  FORMAT  (17F7.3)  16 
READ(1,4)  XHT,CSK,SNK  17 

4  FORMAT  ( 13F9.5)  18 
C   READ  MAP  TITLE»  DATE  OF  ROUTING  COMPUTATION*  MAP  GRID  PLOT  LABEL* 

C        AND  TOTAL  NUMBER  OF  SHIPS  ROUTED 

READ(50*1)  TI »DATE,AL»NST  19 

1  FORMAT  (6A8/8A8,I3)  20 
WRITE(51»2)  NST,KATE  21 

2  FORMAT( 39H1TOTAL  NUMBER  OF  VC2AP3  SHIPS  ROUTED  =  I3/1X3HON  A8*54H  22 
-♦-FROM  12Z  TO  24Z*  AND  ON  FOLLOWING  DAY  FROM  00Z  TC  11Z/)  23 

REWIND  1  24 

DO  80   M=1»NST  25 

IF  (M-l)   10,11,10  26 
C   READ  MAP  GRID  DATA  FOR  DRAW  SUBROUTINE 

10  READ(1»3)  (X(  I  )  , I  =  l»390)» (Y( I )  ,1=1,390)  27 
REWIND  1  28 

C   DRAW  MAP  GRID 

11  CALL  DRAW  (  386 ♦ X , Y » 1 ,0 , LA » I T ♦ 2 . » 2 . »0 ,0 , 2 ♦ 2 *9 , 1 5 ,0 ♦ LAST )  29 
C   READ  SHIP  IDENTIFICATION  NUMBER*  DATE  AND  HOUR  OF  DEPARTURE*  COORDINATES 
C        OF    TRACK  END  POINTS*  CONVERGENCE  FACTOR,  OPTION  TO  PLOT  EARLIER 

C        TRACK*  TIME  STEP  RECIPROCAL,  AND  NUMBER  OF  ITERATIONS 

READ( 50*14)  GL,DATEX,FL,HR,XLG1*XLT1 , XLG2 » XLT2 *RMUL , NN, NSTEP ,LMAX , 30 

+  NP,PALF*PT  31 

14  FORMAT  (A4,2A8*F3.0,F6.1»F5.1,F6.1*F5.1»F3. 0*11, 212. 1 1 , F8 . 5 , F6. 3 )  32 

RSTEP  =  NSTEP  33 
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WRITE(51»15)  LG»KATtX»LF»HR.KATE»XLGl»XLTl»XLG2»XLT2»RMUL 
+  NSTEP,NN,NP 

15  FORMAT ( 15H0ROUTE  OF  SHIP  A4,llH  BEGINS  ON  A8.16H*  JULIAN 
+»1H,/1XF3.0,22H  HOURS  AFTER  1200Z  ON  A8/19H   FROM  LONGITU 
+»16H  AND  LATITUDE  =  F6.1/19H     TO  LONGITUDE  =  F6.1»16H  A 
+DE  =  F6.1//6H  RMUL=F5.0.3X5HLMAX=I2*3X6HNSTEP=I2»3X3HNN=I 
+  1  1//) 

CHECK  ON  OPTION  TO  PLOT  EARLIER  TRACK 
IF  (NN)   16»19»16 

16  READ(50,17)  GL,PL»NK 

17  FORMAT  (2A8.I2) 

READ<50,29)  ( X (  I ) ,  I  =  1 ,20  )  »  ( Y (  I ) » I =1 »20  } 
29  FORMAT  (10F5.2) 

CALL  DRAW  ( NK ♦ X ♦ Y , 2 ♦ 2 » LP » I T , 2 , , 2 . .0 , 0 ♦ 2 , 2 ,9 ♦ 1 5 tO » LAST ) 
WRITE(51»18J  LG,LP 

18  FORMAT(23H0EARLIER  ROUTE  OF  SHIP  A8,15H  ON  JULIAN  DAY  A4/ 
♦BEEN  PLOTTED  USING  PLUS  SIGNS  FOR  SUCCESSIVE  DAILY  POS I T I 

COMPUTATION  OF  GEODESIC  TRACK 

19  ARG  =  (XLG1-10. )/57. 29577951 
C0SLG1=  COSF(ARG) 

SINLG1=  SINF(ARG) 

ARG  =  (XLG2-10. 1/57.29577951 

C0SLG2=  COSF(ARG) 

SINLG2=  SINF(ARG) 

ARG  =  XLT1/57. 29577951 

C0SLT1=  COSF(ARG) 

SINLT1=  SINF(ARG) 

ARG  =  XLT2/57. 29577951 

C0SLT2=  COSF(ARG) 

SINLT2=  SINF(ARG) 

EL  =  SINLT2*C0SLT1*SINLG1 

EM  =-SINLT2*C0SLTl*C0SLGl 


-  C05LT2*SINLT1*SINLG2 
+  C0SLT2*SINLT1*C0SLG2 


EN  =(SINLG2*C0SLG1-C0SLG2*SINLG1 ) *COSLT 1*C0SLT2 
ROOT  =  SQRTF(EL*EL  +  EM*EM  +  EN*EN) 
EL  =  EL/ROOT 
EM  =  EM/ROOT 
EN  =  EN/ROOT 

PR1=  31.205*COSLTl/( 1.+SINLT1 ) 
XI  =  PR1*C0SLG1 
Yl  =  PR1*SINLG1 

PR2=  31.205*COSLT2/( 1.+SINLT2) 
X2  =  PR2*COSLG2 
Y2  =  PR2*SINLG2 
DELX  =  X2  -  XI 
DELY  =  Y2  -  Yl 

S12   =  SQRTF(DELX*DELX  +  DELY*DELY) 
ARC=  S12 

IF  (XLG2-XLG1)  20»21.20 
20  ARG=  ABSFIEN/62.A1  ) 

ARC=  ASINF( ARG*S12 ) /ARG 


»LMAX, 

34 

35 

DATE  =  A436 

DE  =  F6. 

137 

ND  LATITU38 

1,3X3HNP 

=  39 

40 

41 

42 

43 

44 

45 

46 

47 

66H   HAS 

48 

ONS/) 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 

61 

62 

, 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 

81 
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21 


22 
23 

24 
25 

26 
27 


31 

97 
98 


COSA  =  -EN*Y1/31.205  +  EM 

SINA  =   EN*X1/31.205  -  EL 

IF  (COSA)   23.22.23 

ALF  =  SIGNF ( 1.5707963268.SINA) 

GO  TO  27 

ALF  =  ATANF(SINA/COSA) 

IF  (COSA)   24.27.27 

IF  (SINA)   26.25.25 

ALF  =  ALF  +  3.1415926536 

GO  TO  27 

ALF  =  ALF 

N3  =  0 

X(  1  ) 


32 


33 


-  3.1415926536 


Y(l) 
XFIN 
YFIN 
LR  = 
STEP 
TAU 

S 
TVAR 


XI 
Yl 
X2 
Y2 


24. 
16. 
24. 
16. 


YVARS(5 ) 

YVARS(6) 

YVARS( 7) 

Nl 

N2 

DO 

DO 

TC 

DO 


l.O/RSTEP 

0.0 

0.0 

HR/24. 

X(  1  ) 
Y(l  ) 
0.0 


=  1 

=  1 

40   K=2.900 

32   1=1,4 

=  C(  I  )*STEP 

31   J=5,7 


+  TVAR 


YC(J)  =  C( I )*AK( 1-1, J)  +  YVARS(J) 
IF  (ABSF(YC(5)-  9.5)-  7.5)  97,38,38 
IF  (ABSF( YC(6)-16.5)-14.5)  98,38,38 
CALL  TERP 

AP3 

=  ( 16.-YC(6) ) *EN/31.205  +  EM 

=  (YC( 5)-24. )*EN/31.205  -  EL 

=   COSP*CK  +  SINP*SK 


CALL 

COSP 

SINP 

COST 

SINT 

CALL 

DY(5)= 

DY(6)= 


=   SINP*CK  - 

VDERIV 

CAPV*COSP 
CAPV*SINP 


COSP*SK 


DY(7)=  CAPV 

DO  32   J=5,7 

AK( I ,J)  =  STEP*DY( J) 

DO  33   J=5,7 

YVARS(J)  =  (  AK( 1, J)+2.*A<(2,J)+2.*AK( 3,J)+AK(4,J) 

TVAR  =  TVAR  +  STEP 

X(K)  =  YVARS(5) 

Y(K)  =  YVARS(6) 


82 
83 
84 
85 
86 
87 
88 
89 
90 
91 
92 
93 
94 
95 
96 
97 
98 
99 
100 
101 
102 
103 
104 
105 
106 
107 
108 
109 
110 
111 
112 
113 
114 
115 
116 
117 
118 
119 
120 
121 
122 
123 
124 
125 
126 
127 
)/6.  +  YVARS(J)128 
129 
130 
131 


-15- 


Nl   =  K  132 

N2  =  K  133 

IF  (YVARS(7)-ARC)   35.34,34  134 

34  RAT  =  ( ARC-S )/( YVARS(7 )-S)  135 
T  =  STEP*RAT  +  TAU  136 
X(K)  =  (X(K)-X(K-l  )  )*RAT  +  X(K-l)  137 
Y(K)  =  (Y(K  )-Y(K-l  )  )*RAT  +  YU-1)  138 

N2  =  K+l  139 

X(N2)  =  XFIN  140 

Y(N2)  =  YFIN  141 

GO  TO  41  142 

35  S  =  YVARS<7)  143 
TAU=  TAU  +  STEP  144 
IF  (K-900)   36.38,38  145 

36  IF  (ABSF(X(K)-  9.5)-  7.5)  37,38.38  146 

37  IF  (ABSF( Y( Kl-16.5 )-14.5)  40.38,38  147 

38  T  =  TAU  148 
WRITE(51.39)   LG  149 

39  FORMAT(61H0MORE  THAN  899  INTEGRATION  STEPS  OR  WAVE  DATA  FIELD  EXCE150 
+EDED./21H  OTS  ROUTING  OF  SHIP  A4 .4X39HABANDONED  BUT  GEODESIC  TRACK151 
+  IS  PLOTTED/)  152 

N3  =  1  153 

GO  TO  41  154 

40  CONTINUE  155 

41  L  =  0  156 
WRITE(51,42)  157 

42  FORMAT ( 4X 1HL4X2HN16X3HALF7X 1HT7X5HX ( Nl ) 4X4HX F I N5X 5H Y ( Nl ) 4X4HYF I N/ ) 1 58 
PRINT  WEIGHTING  FACTOR  ALPHA  AND  TIME  T  OF  GEODESIC  TRACK 

WRITE(51,43)  L,N1,ALF,T,X(N1 ), XFIN, Y(N1 ), YFIN  159 

43  FORMAT  ( 15, I  6  ,F1 1 . 5 , 5F9  .3  )  160 
ROTATE  AND  TRANSLATE  AXES  TO  PLOT  GEODESIC  TRACK  ON  MAP  GRID 

DO  44   1=1, N2  161 

TEMP  =  .97780241408*X<  I  )  -  . 20952908873* Y ( I )  +  3.0032502718  162 

Y(I)  =  .20952908873*X(  I  )  +  . 97780241408* Y (  I )  -  4.4699538929  163 

44  X( I )  =  TEMP  164 
CALL  DRAW  (  N2 ,X , Y , N3  +  2 ,0 , LG , I T , 2 . , 2 . ,0 , 0 , 2 ,2 , 9 , 1 5 ,0 , LAST )  165 
IF  (N3)   80,45,80  166 

PREPARE  FOR  LMAX  ITERATIONS  TOWARD  MINIMAL-TIME  TRACK 

45  TC  =  HR/24.  167 
X(l)  =  XI  +  24.  168 
Y( 1)  =  Yl  +  16.  169 
YC(5)  =  X( 1  )  170 
YC(6)  =  Y( 1 )  171 
CALL  TERP  172 
DO  9  1=2,399  173 
X( I )  =  0.0  174 

9  Y( I )  =  0.0  175 

X(101)  =  H  176 

YVARS(5)  =  X(l)  177 

YVARS16)  =  Y(l)  178 
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CALL  ANGLE 

Y( 101)  -  XK 

X(20l)  =  XLG1 

Y(201)  =  XLT1 

LR  =  1 

IF  (NP)  81,82,81 

81  ALF  =  PALF 

T   =  PT 
COSA  =  COSF(ALF) 
SINA  =  SINF(ALF) 

82  DO  69   L=1»LMAX 
TVAR  =  HR/24. 
TAU  =  0.0 

Nl  =  XINTFtRSTEP  * 

XN1  =  Nl 

STEP1  =  1.0/RSTEP 

FSTEP  =  -XN1/RSTEP 

Nl  =  Nl  +  2 

DO  46   1=1,10 

46  YVARS( I )  =  0.0 
YVARS(l)  =  1.0 
YVARS(4)  =  1.0 
YVARS(5)  =  X(l) 
YVARS(6)  =  Y(l) 
NK  =  1 

DO  66  K=2,N1 
STEP  =  STEP1 
IF  (K-Nl)   48,47,4 

47  STEP  =  FSTEP 

48  DO  52   1=1,4 

TC  =  C( I )*STEP  +  T 
DO  49   J=l,10 

49  YC( J)  =  C( I )*AK( I- 
IF  (ABSF(YC(5)-  9. 

99  IF  (ABSF(YC(6)-16. 
100  XLAM  =  YC( 1 )*C0SA 
XMU   =  YC(2)*C0SA 
CLAM  =  SQRTF(XLAM* 
CALL  TERP 
CALL  AP3 

ABS  =  (XLAM*CK  +  X 
ORD  =  (XMU  *CK  -XL 
HYP  =  SQRTF(ABS*AB 
VMAJ=  A*ABS/HYP  - 
VMIN=  B*ORD/HYP 

V  =  SQRTF(VMAJ*V 
C0ST=  VMAJ/V 
SINT=  VMIN/V 
C0SP=  CK*C0ST  -  SK 
SINP=  SK*C0ST  +  CK 


T) 


8 


VAR 

It  J)  +  YVARS( J) 
5)-  7.5)  99,65,65 
5)-14.5)100,65,65 
+  YC(3)*SINA 
+  YC(4)*SINA 
XLAM  +  XMU*XMU) 


MU*SK)*A/CLAM 
AM*SK)*B/CLAM 
S  +  0RD*0RD) 
CC 

MAJ  +  VMIN*VMIN) 


*SINT 
*SINT 


179 
180 
181 
182 
183 
184 
185 
186 
187 
188 
189 
190 
191 
192 
193 
194 
195 
196 
197 
198 
199 
200 
201 
202 
203 
204 
205 
206 
207 
208 
209 
210 
211 
212 
213 
214 
215 
216 
217 
218 
219 
220 
221 
222 
223 
224 
225 
226 
227 
228 
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YC(2)*SINP 


YC(4)*SINP 


YC(2 )*YC(3) 


CALL  VDERIV 

FAC1  =  YC(1)*C0SP  + 

DY(1)=  -CAPVX*FAC1 

DY(2)=  -CAPVY*FAC1 

FAC2  =  YC(3)*C0SP  + 

DY(3)=  -CAPVX*FAC2 

DY(4)=  -CAPVY*FAC2 

DY(5)=  CAPV  *  COSP 

DY(6)=  CAPV  *  SINP 

DET  =  YC( 1 )*YC(4)  - 

QUO  =  RBV/CLAM 

FACl=CAPV*QUO*VPBVX/DIV 

FAC2=CAPV*QUO*VPBVY/DIV 

FAC3=  CAPV*QUO*QUO*QUO/DIV 

FAC4=CAPV*DET*QU0/DIV 

D(l)  =  CAPVX»COSP  -  XMU*FAC1 

D(2)  =  CAPVY*COSP  -  XMU*FAC2 

D(3)  =  XMU  *  XMU*FAC3 

D(4)  =  -XLAM*XMU*FAC3 

D(5)  =  -XMU*  DET*FAC3 

DY(7)  =  D(1)*YC(7)  +  D(2)*YC(8) 

D(6)  =  CAPVX*SINP  +  XLAM*FAC1 

D(7)  =  CAPVY*SINP  +  XLAM*FAC2 

D(9)  =  XLAM*XLAM*FAC3 

D(10)=  XLAM*  DET*FAC3 

DY(8)  =  D(6)*YC(7>  +  D(7)*YC(8) 

D(ll)=  (-CAPV*VPBVX*VPBVX/DIV  - 

D(12)=  (-CAPV*VPBVX*VPBVY/DIV  - 


+  D(3)*YC(9)  +  D(4)*YC(10)  +  D(5) 


+  D(A)*YC(9) 
CAPVXX) /QUO 
CAPVXY) /QUO 


+  D(9)*YC( 10) 


-  D( 1 )*YC(9) 
CAPVYY) /QUO 


-  D(2)*YC(9)  -  D(7)*YC(10)  +  D(20) 


)/6.  + 


D(15)=  -VPBVX*FAC4 

DY(9)  =D(11)*YC(7)  +D(12)*YC(8) 

D(17)=  (-CAPV*VPBVY*VPBVY/DIV  - 

D(20)=  -VPBVY*FAC4 

DY( 10)=D( 12)*YC(7)  +D(17)*YC(8) 

DO  52   J=l,10 

52  AK( I »J )  =  STEP  *  DY( J) 
DO  53   J=l,10 

53  YVARS(J)  =  (  AK(l.J)+2«*AK(2tJ)+2.*AK(3»J)+AK(4.J) 
TVAR  =  TVAR  +  STEP 

TAU   =  TAU  +  STEP 
IF  (Nl-K)   54.56.54 

54  IF  (LMAX-L)   62.55.62 

55  IF  (  (K-D/NSTEP  +  1-NK)   62.62.56 

56  NK  =  NK  +  1 
YC(5)=  YVARS(5) 
YC(6)=  YVARS(6) 
LR  =  0 

CALL  TERP 

CALL  AP3 

LR  =  1 

IF  (LMAX-L)   61.60,61 


-  D(6)*YC( 10) 


229 
230 
231 
232 
233 
234 
235 
236 
237 
238 
239 
240 
241 
242 
243 
244 
245 
246 
247 
248 
249 
250 
251 
252 
253 

+  D( 10)254 
255 
256 
257 

+  D( 15)258 
259 
260 
261 
262 
263 
264 
YVARS( J)265 
266 
267 
268 
269 
270 
271 
272 
273 
274 
275 
276 
277 
278 
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60  X(NK)  =  YVARS(5) 
Y(NK)  =  YVARS(6) 
X(NK+100)  =  H 
XINK+300)  =  TAU 
CALL  ANGLE 
Y(NK+100)  =  XK 
X(NK+200)  =  XLG 
Y(NK+200)  =  XLT 

61  IF  (Nl-K)   62*67*62 

62  DELX  =  YVARS(5)  -  X(  1  ) 
DELY  =  YVARS(6)  -  Y(  1) 
IF  (DELX*DELX  +  DELY*DELY  -  S12*S12)  63.65*65 

63  IF  ( ABSF( YVARS( 5)-  9.5)-  7,5)  64,65,65 

64  IF  (ABSF(YVARS(6)-16.5)-14.5>  66*65*65 

65  Nl  =  < 
T   =  TAU 
GO  TO  56 

66  CONTINUE 
PRINT  ALPHA*  T,  X*  AND  Y  AT  END  OF  EACH  ITERATION 

67  WRITE(51»43)  L ♦ Nl » AL F * T ♦ Y VARS ( 5 ) »XF I N » YVARS ( 6 ) » YF I N 
XLAM  =  YVARS(l)*COSA  +  YVARS ( 3 ) *S I NA 

XMU  =  YVARS(2) *COSA  +  YVARS ( 4 ) *S I NA 
CLAM  =  SQRTF(XLAM*XLAM  +  XMU*XMU) 
ABS  =  (XLAM*CK  +   XMU*SK) *A/CLAM 
ORD  =  (  XMU*CK  -  XLAM*SK)*B/CLAM 
HYP  =  SQRTF(ABS*ABS  +  ORD*ORD) 
VMAJ=  A  *  ABS/HYP  -  CC 
VMIN  =  B  *  ORD/HYP 
DELX  =  YVARS(5)  -  24. 
DELY  =  YVARS(6)  -  16. 

EMFI  =  (DELX*DELX  +  DELY*DELY  +  973. 75 )/ 1043 .638743 
XDOT=  (CK*VMAJ-SK*VMIN)*EMFI/ 8. 566041666  7 
YDOT=  (SK*VMAJ+CK*VMIN)*EMFI/8.566041666  7 
DIFX  =  XFIN  -  YVARS( 5  ) 
DIFY  =  YFIN  -  YVARS(6) 
DET  =  XDOT*YVARS(8)  -  YDOT*YVARS ( 7 ) 
DIFT=  (YVARS(8)*DIFX  -  YVARS ( 7 ) *DI FY ) /DET 
DIFA=  (XDOT*DIFY  -  YDOT*D I FX ) / DET 

T    =  DIFT  +  T 
ALF   =  DIFA/RMUL  +  ALF 
COSA  =  COSF(ALF) 
SINA  =  SINF(ALF) 
PRINT  NEW  VALUES  OF  ALPHA  AND  T 

WRITE(51»50)  ALF»T  320 

50  FORMAT  ( 1 1XF 1 1 . 5 » F9. 3  )  321 

69  CONTINUE  322 
IF  (DIFX*DIFX  +  DIFY*DIFY  -  EMF I *EMF I  * . 2366 )   72,72,70  323 

70  WRITE(51,71  )  LG  324 

71  FORMAT(20H0  OTS  ROUTE  OF  SHIP  A4,47H  MORE  THAN  100  MILES  FROM  DEST325 
+  INATION  BUT  TRACK/69H   IS  PLOTTED.  INCREASE  RMUL  OR  LMAX,  OR  BOTH, 326 
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+  TO  IMPROVE  CONVERGENCE.)  327 
TABOLATE  FINAL.  TRACK  DAILY  POSITION*  WAVE  HEIGHT  AND  DIRECTION 

72  WRITE(51»73)  328 

73  FORMAT  (  1  HO /^X4HDA YS 7X 5HLONG I 4X 5HL AT  I -5X4HWA VE5X14HW AVE  DIRECTION/329 
+2X9HOF  TRAVEL4X5H-TUDE4X4HTUDE5X6HHEIGHT6X10HFROM  NORTH/)  330 

WRITE(51»74)<X(K+3  00),X(K+200),Y(K+200),X(K+100),Y(K+100),K=1,NK)  331 

74  FORMAT  (F9.2,F11. 1 ,F8.1,P9.0,F14.0)  332 
ROTATE  AND  TRANSLATE  AXES  FOR  PLOT  OF  DAILY  POSITIONS 

DO  7  5   1=1, NK  333 

TEMP  =  .97780241408*X( I )  -  . 20952908873* Y ( I )  +  3.0032502718  334 

Y(I)  =  .20952908873*X( I )  +  . 97780241408*Y ( I )  -  4.4699538929  335 

75  X( I )  =  TEMP  336 
PUNCH  11  CARDS  USEABLE  FOR  A  LATER  PLOT  OF  TRACK 

WRITE(52,76)  LG,LF,NK  337 

76  FORMAT  ( 2 A8 , I  2 , 6 1  X  1H  )  338 
WRITE(52,78)  (  (  (  RX ( I , J ) , I  =  1 , 1  0  )  » LG ♦ LF , J  )  ♦ J  =  1  ,  2  )  ,  339 

+  (  (  (RY(  I  ,J  )  ,  1  =  1,10)  ,LF,LG, J)  ,J  =  1  ,2  )  340 

78  FORMAT  ( 1  OF  5. 2 , 2A8 , 1  2 , 1 IX 1H  )  341 

CALL  DRAW  ( NK , X ♦ Y , 3 ♦ 4 , LF , I T , 2 . , 2 . ,0 ,0 » 2 , 2 ,9 , 1 5 ,0 , LAST )  342 

WRITE(51»93)  343 

93  FORMAT  (1H1)  344 

PROCEED  TO  COMPUTE  THE  ROUTE  OF  NEXT  SHIP 

80  CONTINUE  345 

STOP  346 

END  347 
SUBROUTINE  TERP 

DIMENSION  HT(4,4)»CT(4,4),ST(4,4)»P{4),Q(4),PX(4),QY(4),PXX(4),QYY1 

+  (4)  ,HD(4)  ,CD(4) ,SD(4)  ,HS(4) ,CS(4) ,SS<4)  ,HP(4) ,CP(4)  ,SP(4)  ,HXS(4),  2 

+  CXS(4)  ,SXS(  4 ) ,HPX(4)  ,HPY(4) , HPXX ( 4 ) ,HPXY ( 4 )  ,HPYY(4) ,CPX(4) ,CPY(4)  ,3 

+  CPXX(4)  ,CPXY(4)  ,CPYY(4)  ,SPX(4)  ,SPY(4)  , SPXX ( 4  )  , SPXY ( 4  )  »SP YY ( 4 ) ,C ( 4 ) 4 

COMMON     YC( 10) ,LR,A,B,CC,H»CK»SK  5 

COMMON/L1/XHT(5  760) ,CSK(5760) ,SNK( 5760) ,TC  6 

COMMON/L2/HX»HY,HXX,HXY,HYY  7 

COMMON/ L3/DKX,DKY»DKXX,DKXY,DKYY  8 

DTC  =  2.*TC  9 

L  =  XINTF(DTC)  10 

IF  (L-3)   1,1,7  11 

1  TT  =  (-INTF(DTC)+DTC)»2.  -  1.  12 
TP1=  TT  +  1.  13 
TM1=  TT  -  1.  1* 
T2M=  TP1*TM1  15 
IF  (L)   2,2,3  16 

2  K4  =  3  17 
TM3=  TT  -  3.  18 
C(l)=  TMl*TM3/8.  19 
C(2)=-TPl*TM3/4.  20 
C(3)=  T2M/8.  21 
GO  TO  16  22 

3  K4  =  4  23 
IF  (L-2)   4,4,6  24 
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4  G  =  ( 3.*TT+2. )*TT  -  9. 
F  =  -4.*TT  +  G 

C(  1  )  =  -T2M*TM1/16. 
C(2)=  G*TM1/16. 

5  C(3)=  -F*TP1/16. 
C(4)  =  T2M*TP1/16. 
GO  TO  15 

6  C(l)=  -T2M*TM1/16. 
C(2)=((2.*TT+1.)*TT-7.)*TM1/12. 
C(3)=(  ( l.-TT)*TT+4.)*TPl/8. 
C<4)  =  T2M*TPl/48. 

GO  TO  15 

7  L  =  XINTFCTC-2.)  +  4 
IF  (L-10)   9.8,8 

8  <4  =  1 

L  =  9 
C(l)=  1. 
GO  TO  16 

9  TT  =  (-INTF(TC)+TC)*2.  -  1. 
TP1=  TT  +  1. 

TM1=  TT  -  1. 
T2M=  TP1*TM1 

G  =  ( 3.*TT+2. )*TT  -  9, 

F  =  -4.*TT  +  G 
C( 1 )=  -T2M*TM1/16. 
IF  (L-9)   11,10,8 

10  ^  =  2 

C(2)=  (5*TM1  +  (T2M-F)*TP1 )/l6. 
GO  TO  15 

11  C(2)=  G*TM1/16. 

IF  IL-8)   13. 12,10 

12  K4  =  3 

C(3)=  ( T2M-F)*TP1/16. 
GO  TO  15 

13  K4  =  4 

IF  (L-5)   14,5,5 

14  C(l)=  -T2M*TMl/6. 
C(2)=((5.*TT+2.)*TT-11.)*TM1/16. 
C(3)=((-5.*TT+4.)*TT+13.)*TPl/24. 
C(4)  =  T2M*TP1/16. 

L-l 

XINTF( YC(5) )  -  2 

XINTF( YC(6) )  -  2 

(-INT- (YC(5)  )+YC( 5)  )*2.0  -  1. 

(-INTF( YC(6)  )+YC(6)  )*2.0  -  1. 

XX  +  1.0 

XX  -  1.0 

YY  +  1.0 

YY  -  1.0 

XP1*XM1 


15 

L  = 

16 

M  = 

M  = 

XX  = 

YY  = 

XP1  = 

XM1  = 

YP1  = 

YM1  = 

X2M^ 

25 
26 
27 
28 
29 
30 
31 
32 
33 
34 
35 
36 
37 
38 
39 
40 
41 
42 
43 
44 
45 
46 
47 
48 
49 
50 
51 
52 
53 
54 
55 
56 
57 
58 
59 
60 
61 
62 
63 
64 
65 
66 
67 
68 
69 
70 
71 
72 
73 
74 
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17 


18 


19 


20 


Y2M=  YP 
P(l)  = 
P(2)  = 
P(3)  = 
P  (  4  )    = 
kll)    = 
6(2)  = 
0(  3)  = 
8(4 )  - 
PX(4)= 
PX(  1  )  = 
PX(2>=< 
PX(3)= 
&Y<4)= 
OY (  1  )  = 
&Y(2>=( 
JY( 3)= 
IF  (LR) 
PXX(4) = 
^XX( 1)= 
PXX( 2) = 
PXX(3)= 
wYY<4)  = 
JYY (  1  )  = 
UYY<2)= 
CJYY  (3)  = 
DO  2  7 
HP  (  K  )  = 
CP(K)  = 
5P(K)  = 
HPX(K) = 
HPY (K) = 
CPX(K) = 
CPY(K)= 
SPX(K) = 
SPY (K)  = 
IF  (LR) 
HPXX(<) 
HPXY(K) 
HPYY(K) 
CPXX( K) 
CPXY(l^  ) 
CPYY(K) 
SPXX<<) 
SPXYU ) 
SPYY(K) 
KK  =  ( 
DO  23 
HD(J)  = 
CD(J)  = 


1*YM1 

(XX+2. )*X2M*XMl*XMl/32. 

( ( ( ( -3.*XX*XX+12. )*XX-2. )*XX-2  5.)*XX+18.)/32. 

(-XX*XX+9. ) /8.  -  P(2) 

(-XX+2. )*X2M*XPl*XPl/32. 

(YY+2. )*Y2M*YMl*YMl/32. 

( ( ( (-3.*YY*YY+12.)*YY~2. )*YY-25. )*YY+18. )/32. 

(~YY*YY+9. )  /8.  -  Q( 2 ) 

(-YY+2. )*Y2M*YPl*YPl/32. 

( (-5.*XX+10.) *XX-3. )*XP1*XP1/16. 

XX/2.  -  PX(4) 
( {-15.*XX*XX+36. )*XX-4. )*XX  2  5.  )  /16. 
-XX/2.  -  PX(2) 
(  t-5.*YY+10,  )  *YY-3.  )*YP1*YP1/16. 

YY/2.  -  QY(4) 
l(-15.*YY*YY+36.  )*YY-4,  )  *YY  2  5.  ) /16. 
-YY/2.  -  OY(2  ) 

17»18»17 
(  (-5.*XX*XX  +  6.  ) *XX+1. ) /2. 

1.  -  PXX(4  ) 
((-15.*XX*XX+18. )*XX-1. ) /2. 

-1.  -  PXX( 2) 
( (-5.*YY*YY+6. )*YY+1. ) /2. 

1.  -  QYY(4) 
(<-15.*YY*YY+18. )*YY-1. )/2« 

-1.  -  QYY(2 ) 
K=1»K4 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

19,20,19 

=  0.0 

=  0.0 

=  0.0 

=  0.0 

=  0.0 

=  0.0 

=  0.0 

=  0.0 

=  0.0 
(K+L )*32+N) *18  +  M  -  594 
J=l,4 

0.0 

0.0 


75 

76 

77 

78 

79 

80 

81 

82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 

100 

101 

102 

103 

104 

105 

106 

107 

108 

109 

110 

111 

112 

113 

114 

115 

116 

117 

118 

119 

120 

121 

122 

123 

124 
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21 


11 


23 


24 


25 


5D(J)  = 
HS(J)  = 
CS(J)  = 
c>S(J)  = 
IF  (LR) 
HXS( J) 
CXS( J) 


26 


0.0 

0.0 

0.0 

0.0 

21.22.21 
=  0.0 
=  0.0 


SXS(J)  =  0.0 
JJ  =  J*18  + 

DO  2  3 
II  =  I 

HT( I »J) 

CT(  I  »J) 

ST(  I  »J) 

DO  25 


11 


K.K 
1  =  1.4 
+  JJ 

=  XHT(  II  ) 
=  CSK( II) 
=  SNK< I  I ) 
1  =  1.4 
DO  25   J=l»4 
HD( I )  =  Q( J)*HT( I » J) 
CD( I )  =  Q( J)*CT( I » J) 
i>D(  I  )  =  Q( J)*ST< I.J) 
HS(  I )  =  P( J)*HT< J.  I  ) 
CS(  I  )  =  P( J)*CT(J.I  ) 
SS(  I  )  =  P( J)*ST(J.  I  ) 
IF  (LR)  24.25.24 
HXS( I )=PX( J)*HT(J» I ) 
CXS( I )=PX( J)*CT(J» I  ) 
^XS( I )=PX( J)*ST( J. I ) 
CONTINUE 
DO  11       1=1,4 
HP(K)  =  HD( I )*P( I )  + 
CP(K)  =  CD( I )*P( I )  + 
SP(K)  =  SD( I)*P( I )  + 
HPX(K)=HD( I )*PX( I )  + 
CPX(K)=CD( I )*PX( I )  + 
SPX(K)=SD( I )*PX( I )  + 
HPY(K)=HS( I )*QY( I )  + 
CPY(K)=CS( I )*QY( I )  + 
bPY(K)=SS( I )*QY( I )  + 
IF  (LR)26, 27.26 
HPXX(K)=  HD( I )*PXX( I 
HPXY(<)=HXS( I )*  QY( I 
HPYY(K)=  HS( I )*QYY( I 
CPXX(K)=  CD( I )*PXX( I 
CPXY(K)=CXS( I )*  QY( I 
CPYY(K)=  CS( I )*QYY( I 
SPXX(K)=  SD( I )*PXX( I 
SPXY(K) =SXS( I )*  QY( I 
SPYY(K)=  SS( I )*QYY( I 
CONTINUE 
H   =  0.0 
CK  =  0.0 


HD(  I  ) 
CD(  I  ) 
SD(  I) 
HS(  I  ) 
CS(  I  ) 
SS(  I  ) 


+  HXS( I ) 
+  CXS( I ) 
+  SXS( I  ) 


HP(K) 

CP(K) 

SP(K) 

HPX(K) 

CPX(K) 

SPX(K) 

HPY(K) 

CPY(IC) 

SPY(K) 


HPXX(K) 
HPXY(K) 
HPYY(K) 
CPXX(K) 
CPXY(K) 
CPYY(K) 
SPXX(IC) 
SPXY(K) 
SPYY(K) 


125 
126 
127 
128 
129 
130 
131 
132 
133 
134 
135 
136 
137 
138 
139 
140 
141 
142 
143 
144 
145 
146 
147 
148 
149 
150 
151 
152 
153 
154 
155 
156 
157 
158 
159 
160 
161 
162 
163 
164 
165 
166 
167 
168 
169 
170 
171 
172 
173 
174 


-23- 


SK  =  0,0 

HX  =  0.0 

HY  =  0.0 

CKX  =  0.0 
CKY  =0.0 

SKX  =  0.0 

5KY  =  0.0 
IF  (LR)   28*29.28 

28  HXX  =  0.0 
HXY  =  0.0 
HYY  =  0.0 
CKXX=  0.0 
CKXY  =  0.0 
CKYY=  0.0 
SKXX=  0.0 
SKXY=  0.0 
SKYY  =  0.0 

29  DO  31  K=1,K4 

H   =  C(K)*HP( K  )  + 

CK   =  C(K)*CP( K)  + 

SK   =  C(K)*SP(K)  + 

HX  =  C (K)*HPX( K) 

HY  =  C(K)*HPY(K) 

CKX  =  C(K)*CPX(K) 

CXY  =  C(K)*CPY(K) 

SKX  =  C(K)*SPX(K) 

S<Y  =  C(K)*SPY(K) 
IF  (LR)   30,31,30 

30  HXX  =  C(K)*HPXX(K) 
HXY  =  C(K)*HPXY(K) 
HYY  =  C(K)*HPYY(K) 
CKXX  =  C(K)*CPXX(K) 
CKXY  =  C(K)*CPXY(K) 
CKYY  =  C(K)*CPYY(K) 
SKXX=  C<K)*SPXX(K) 
S,<XY  =  C(K>*SPXY(K) 
S.<YY  =  C(K)*SPYY(K) 

31  CONTINUE 


32 


33 


H 

CK 

SK 


HX 

HY 

CKX 

CKY 

SKX 

SKY 

HXX 

HXY 

HYY 

CKXX 

CKXY 

CKYY 

SKXX 

SKXY 

SKYY 


RAD  =  SQRTF(CK*CK  +  SK*SK) 

CK   =  CK/RAD 

SK   =  SK/RAD 

DKX  =  CK*SKX  -  SK*CKX 

DKY  =  CK*SKY  -  SK*CKY 

IF  (LR)   32,33,32 

DKXX=  CK*SKXX-SK*CKXX 

DKYY=  CK*SKYY-SK*CKYY 

DKXY=  CK*SKXY-SK*CKXY 

RETURN 

END 


+  CKY*SKX  -  SKY*CKX 


175 
176 
177 
178 
179 
180 
181 
182 
183 
184 
185 
186 
187 
188 
189 
190 
191 
192 
193 
194 
195 
196 
197 
198 
199 
200 
201 
202 
203 
204 
205 
206 
207 
208 
209 
210 
211 
212 
213 
214 
215 
216 
217 
218 
219 
220 
221 
222 
223 
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SUBROUTINE  ANGLE 

COMMON     YC(  10) .LR.A.B.CCH.CK.SK  1 

COMMON/L6/YVARS( 10) .XK.XLG.XLT  2 

DELX  =  YVARS(5)  -  24.  3 

DELY  =  YVARS(6)  -  16.  4 

COSXK=  -DELX*CK  -  DELY*SK  5 

5INXK=   DELX*SK  -  DELY*CK  6 

IF  (COSXK)   2.1»2  7 

1  XK  =  SIGNF(90..SINXK)  8 
GO  TO  6  9 

2  XK  =  ATANF(SINXK/COSXK)*57. 29577951  10 
IF  (COSXK)   3.6.6  11 

3  IF  (SINXK)   5*4.4  12 

4  XK  =  XK  +  180.  13 
GO  TO  6  14 

b  XK  =  XK  -  180.  15 

to  IF  (XK)   7.8.8  16 

7  XK  =  360.  +  XK  17 

8  XT  =  DELX*. 98480775  -  DEL Y*. 173648 18  18 
YT  =  DELX*. 17364818  +  DELY*. 98480775  19 
RAD=  SQRTF(XT*XT  +  YT*YT)  20 
IF  (XT)   10.9.10  21 

9  XLG  =  SIGNF(90.0.YT  )  22 
GO  TO  14  23 

10  XLG  =  ATANF(YT/XT)*57. 29577951  24 
IF  (XT)   11.14,14  25 

11  IF  (YT)   13.12.12  26 

12  XLG  =  XLG  +  180.  27 
GO  TO  14  28 

13  XLG  =  XLG  -  180.  29 

14  XLT  =  -ATANF(RAD/31.205)*114. 591559  +  90.0  30 
RETURN  31 
END  32 
SUBROUTINE  VDERIV 

COMMON     YC(10)  .LR.A.B.CCH.CK.SK.  1 

COMMON/ L3/ DKX.DKY.DK XX .DKXY.DKYY  2 
C0MM0N/L4/AX.AY.BX.BY.CX.CY.AXX.AXY.AYY.BXX.BXY.BYY.CXX.CXY.CYY    3 

COMMON/ L5 /COST. SI  NT »V .CAP V. CAP VX. CAP VY.CAPVXX.CAPVXY.CAPVYY.  4 

+  VPBVX»VPBVY,DIV.RBV,CAPVP  5 

SIN2T  =  2.*SINT*COST  6 

COS2T  =  2.*COST*COST  -  1.  7 

DELX  =  YC(5)  -  24.  8 

DELY  =  YC(6)  -  16.  9 

EMFI  =  (DELX*DELX  +  DELY*DELY  +  973. 75  J / 1043 . 638743  10 

EMFIX=  DELX/521. 8193715  11 

EMFIY=  DELY/521. 8193715  12 

AMCX  =  A*AX  -  CC*CX  13 

AMCY  =  A*AY  -  CC*CY  14 

AMC2  =  A*A  -  CC*CC  15 
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BCA2  =  B*B  -  AMC2  16 

BCASN2=  BCA2*SIN2T  17 

BCACS2=  BCA2*COS2T  18 

ROOT  =  SQRTF((B*B  +  AMC2  +  BCACS2)/2.)  19 

BCACOS=  B*CC*COST/A  20 

BCASIN=  B*CC*SINT/A  21 

REC  =  BCACOS  +  ROOT  22 

FP  =  .5*BCASN2/ROOT  +  BCASIN  23 

VPBV   =  FP/REC  24 

IF  (LR)  2»1.2  25 

1  V   =  (B*AMC2/REC) /A  26 

2  VSR  =  V/8. 5660416667  27 
CAPV  =  VBR*EMFI  28 

FX  =  (CC*BX  +  B*CX  -  B*CC*AX/A)/A  29 

FY  =  (CC*BY  +  6*CY  -  B*CC*AY/A)/A  30 

RX  =  AMCX*SINT*SINT  +  B*BX*COST*COST  +  . 5*BCASN2*DKX  31 

RY  =  AMCY*SINT*SINT  +  B*BY*CCST* COST  +  •  5*BC ASN2*DKY  32 

RECX=  FX*COST  +  BCASIN*DKX  +  RX/ROOT  33 

RECY=  FY*COST  +  BCASIN*DKY  +  RY/ROOT  34 

FACX=  2.*AMCX/A;v,C2  +  BX/B  -  AX/A  -  RECX/REC  35 

FACY=  2.*AMCY/AMC2  +  BY/B  -  AY/A  -  RECY/REC  36 

CAPVX  =  (FACX*EMFI  +  EMFIX)*V3R  37 

CAPVY  =  (FACY*EMFI  +  EMFIY)*V3R  38 

IF  (LR)   4»3.4  39 

3  CAPVP  =  VPBV  *  CAPV  40 
RETURN  41 

4  AMCXY  =  AX*AY  -  CX*CY  +  A*AXY  -  CC*CXY  42 
AMCXX  =  AX*AX  -  CX*CX  +  A*AXX  -  CC*CXX  43 
AMCYY  =  AY*AY  -  CY*CY  +  A*AYY  -  CC*CYY  44 
FXY=(CC*BXY  +  B*CXY  -  B*CC*AXY/A  +  BX*CY  +  CX*BY- ( B*CY+CC*BY ) *AX/A45 

t     -CC*BX*AY/A  -  B*CX*AY/A  +  { 2 .*B*CC*AX* AY / A ) /A ) / A  46 

FXX=(CC*BXX  +  B*CXX  -  B*CC*AXX/A  +  2.*BX*CX  -  ( B*CX+CC*BX ) *2 .*AX/ A47 

+     +( 2.*B*CC*AX*AX/A) /A) /A  48 

FYY=(CC*BYY  +  B*CYY  -  B*CC*AYY/A  +  2.*BY*CY  -  ( B*CY+CC*BY ) *2 .*A Y/A49 

+     +( 2.*B*CC*AY*AY/A) /A) /A  50 

BCAX  =  8*BX  -  AMCX  51 

BCAY  =  B*BY  -  AMCY  52 

RXY  =  AMCXY*SINT*SINT  +  ( BCAX*DKY+BCAY*DKX ) *S I N2T  -  BCACS2*DKX*DKY53 

+        (BX*BY+B*BXY)*COST*COST  +  .  5*BCASN2*DKXY  54 

RXX  =  AMCXX*SINT*S1NT  +  ( 2 «*BCAX*S I N2T-BCACS2*DKX ) *DKX  55 

+        (6X*BX+B*BXX)*C0ST*C0ST  +  . 5*BCASN2*DKXX  56 

RYY  =  AM.CYY*SINT*SINT  +  (  2  .*BC  AY*S  I  N2T-BCACS2*DKY  )  *DKY  57 

+        (BY*BY+B*BYY)*COST*COST  +  .  5*BCASN2*DKYY  58 

RECXY  =  FXY*COST  +  ( FX*DKY  +  FY*DKX ) *S I  NT  -  BCACOS*DKX*DKY  59 

+          BCASIN*DKXY  +  ( (-RX*RY/ROOT)/ROOT+RXY)/ROOT  60 

RECXX  =  FXX*COST  +  ( 2 .*FX*S INT-BCACOS*DKX ) *DKX  61 

+          BCASIN*D£XX  +  { (-RX*RX/ROOT )/ROOT+RXX) /ROOT  62 

KECYY  =  FYY*COST  +  (  2  .  *FY*S  I  NT-BCACOS*D(CY  )  *DKY  63 

+          BCASIN*DKYY  +  (  (  -RY*R Y/ROOT  )  /ROOT+RYY  )  /ROOT  64 

FACXY  =  (-2.*AMCX*AMCY/AMC2+AMCXY)*2./AMC2  +  ( -BX*BY/B+BXY ) /B  65 
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f  (AX*AY/A-AXY)/A  +  ( RECX*REC Y/R EC-RECXY ) / REC  66 

FACXX  =  (-2.*AMCX*AMCX/AMC2+AMCXX)*2./AMC2  +  ( -BX*BX/B+BXX ) /B  67 

*-  (AX*AX/A-AXX)/A  +  ( RECX*RECX/REC-RECXX  )  / REC  68 

FACYY  =  (-2.*AMCY*AMCyVAMC2+AMCYY)*2./AMC2  +  ( -BY*B Y/B+B YY ) /B  69 

+  (AY*AY/A-AYY)/A  +  ( REC Y*RECY/REC-RECYY ) / REC  70 

CAPVXY  =  ( (FACX*FACY+FACXY)*EMFI  +  FACX*EMFIY  +  FAC Y*EMF I  X ) *VBR  71 
CAPVXX  -  ( (FACX*FACX+FACXX )*EMFI  +2 • *FACX*EMF I X  +  . 00 1916371938 ) *VBR72 
CAPVYY  -  ( ( FACY*FACY+FACYY)*EMFI  +2 . *FAC Y*EMF I Y+. 00 1 91 637 1938 ) * VBR73 
r'PX  =  (  i-.5*6CASN2*RX/ROOT  )/ROOT  +  BCAX*SIN2T  -  BCACS2*DKX ) /ROOT   74 

+        FX*SINT  -  BCACOS*DKX  75 
FPY  =  ( (-.5*BCA5N2*RY/ROOT )/ROOT  +  6CAY*SIN2T  -  BCACS2*DKY ) /ROOT   76 

+        FY*5INT  --  BCACOS*DKY  77 

VPBVX  =  (-FP*RECX/REC  +  FPXJ/REC  78 

VPBVY  =  (-FP*RECY/REC  +  FPY) /RFC  79 

RECP=  -.5*BCASN2/ROOT  -  BCASIN  80 

FPP  =  ( ( ,25*BCASN2*BCASN2/ROOT)/ROOT  +  BCACS2)/ROOT  +  BCACOS  81 

VP8VP  =  (-FP*RECP/REC  +  FPP) /REC  82 

DIV  =  VP6V*VPBV  -  VP3VP  +  1.  83 

RBV  =  SGRTF <VPBV*VP6V  +1.)  84 

RETURN  85 

END  86 
SUBROUTINE  AP3 

COMMON     YC( 10) ,LR»A,B»CC,H,CK»SK  1 

COMMON/L2/HX,HY»HXX»HXY»HYY  2 

C0MM0N/L4/AX,AY,BX,BY,CX»CY,AXX,AXY,AYY,BXX,BXY,BYY,CXX,CXY,CYY  3 

Rl  =  SQRTF(  (  .062760850324*H-.60018313990)*H+4. 7014047597)  A 

VF  =  0.021541997619*H  +  19.278272298  -  Rl  5 

R2  =  SQRTF(  (  .060104035000*H-.96636105838)*H+6. 1294779871 )  6 

B  =  -0.12663045716*H  +  19.585778258  -  R2  7 

DVFM=  (-.  062760850324*H+.  30009156995  )/Rl  8 

DVF  =  DVFM  +  .021541997619  9 

DbK    =  (-.060104035000*H+.483l8052919)/R2  10 

DS  =  DBM  -  0.12663045716  11 

D2VF  =  (DVFM*DVFM  -  . 062760850324 ) /Rl  12 

D2B  =     (DBM  *  DBM  -  . 060104035000 ) /R2  13 

IF  (H-17.)   1,1,2  14 

.  R3  =SQRTF( ( ,083601632403*H-1. 3340008783 )*H+7. 1705253492)  15 

VH  =  -0.24791650490*H  +  19.793624009  -  R3  16 

DVHM=  (-.083601632403*H+.66700043915)/R3  17 

DVH  =  DVHM  -  0.24791650490  18 

D2VH=  (DVHM*DVHM  -  . 08 360 163240 3 ) /R3  19 

50  TO  3  20 

I    R4  =SQRTF( ( .055777533214*H-3.0851911409)*H+45. 698170763)  21 

VH  =  -0.31013284648*H  +  14.848653764  +  R4  22 

DVHM=  (  .055777533214*H-1.54259557045)/R4  23 

DVH  =  DVHM  -  0.31013284648  24 

D2VH=  (-DVHM*DVHM+  . 05 57775332 14 ) /R4  25 

J   A  =  (VF+VH)*.5  26 

CC  =  A-VH  27 

DA  = (DVF+DVH)*.5  28 
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DC  =DA-DVH  29 

AX  =  DA*HX  30 

AY  =  DA*HY  31 

BX  =  DB*HX  32 

BY  =  DB*HY  33 

CX  =  DC*HX  34 

CY  =  DC*HY  35 

IF  (LR)   4»8t4  36 

4  XX  =  HX*HX  37 

XY  =  HX*HY  38 

YY  =  HY*HY  39 

D2A  =( D2VF+D2VH)*. 5  40 

D2C  =  D2A  -  D2VH  41 

AXX=  DA*HXX  +  D2A*XX  42 

AXY=  DA*HXY  +  D2A*XY  43 

AYY=  DA*HYY  +  D2A*YY  44 

bXX=  DB*HXX  +  D2B*XX  45 

BXY=  DB*HXY  +  D2B*XY  46 

BYY=  DB*HYY  +  D2B*YY  47 

CXX=  DC*HXX  +  D2C*XX  48 

CXY=  DC*HXY  +  D2C*XY  49 

CYY=  DC*HYY  +  D2C*YY  50 

8  RETURN  51 

END  52 
END 
FINIS 

EXECUTER. 

OB  0574         BLEICK  BOX  6 

C2AP3           APRIL  14  1967              JUL26,66            1 

666JUL26»66J207     00.  154.0  41.0-123.0  38.0   10  610 

COOP,  BOX  6,  BLEICK  , I / 1/0/2/S/56/ 57 » 5 , 10000  ,  0  ,  TAPE  -  14  APR  67. 

BINARY, 56. 

RELOCOM. 

FTN,L,A,E. 

PROGRAM  TAPE 

DIMENSION  ND(3969)  1 

COMMON     X(390 ) ,Y(390 ) ,MD(63,63) ,  2 

+           XHT(18»32, 10) ,CSK( 18,32,10) ,SNK(18, 32, 10)  3 

EQUIVALENCE  ( MD ,ND ) , ( I D ,ARG > , ( I h,H )  4 

REWIND  1  5 

REWIND  2  6 

ISCALE  =  2000000000000000B  7 

READ  COORDINATES  FOR  MAP  GRID  OF  DRAW  SUBROUTINE 

READ(5U,1)  (X(  I  ), 1  =  1,390) ,  ( Y (  I )  ,  I  =  1  ,  390 )  8 

1  FORMAT  (15F5.3)  9 

WRITE (2 ,7)  X,Y  10 

7  FORMAT  (17F7.3)  11 

IF  (IOCHECK.2)  2,4  12 
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2  WRITE(  51»3)  13 

3  FORMAT  (37HU   PARITY  ERROR  OCCURRED  ON  X»Y  WRITE/)  14 

4  DO  43   K=ltlO  15 
READ  WAVE  DIRECTION  FROM  FLEET  NUMERICAL  WEATHER  FACILITY  TAPE 

BUFFER  IN(1,2)  ( ND ( 1 ) »ND( 3969 ) )  16 

b    IF(UNIT.l)  5.14.8.10  17 

8  WRITE( 51.9)  K  18 

9  FORMAT  (44H0   DIRECTION  EOF  OR  EOT  ERROR  OCCURRED  ON  K=I3/)  19 
GO  TO  6  20 

10  WRITE(51.11)  K.  21 

11  FORMAT  (49HG  DIRECTION  PARITY  OR  LENGTH  ERROR  OCCURRED  ON  K«I3/)  22 
M  =  LENGTHF(l)  2  3 
IF  (M-3969)  12.14,12  24 

12  WRITE ( 51,13)  M  25 

13  FORMAT  (28H0   DIRECTION  BUFFER  LENGTH  =16/)  26 
COMPUTE  COSINE  AND  SINE  OF  WAVE  DIRECTION  K  MEASURED  FROM  X  AXIS 

14  DO  17  1=1.18  27 
DELX  =  1-24  28 
DO  17  J=l,32  29 
DELY  =  J-16  30 
ROOT  =  SQRTF(DELX*DELX  +  DELY*DELY )  31 

ID  =  MD( I+8.J+16) /2046  +  ISCALE  32 

ARG  =  (ARG  +  0.0)*11. 17010721  33 

COS  =  COSF(ARG)  34 

SIN  =  SINF(ARG)  35 

CSMI.J.K)  =(-DELX*COS  -  DELY*SIN ) /ROOT  36 

17  SNK(I,J,K)  =  (DELX*SIN  -  DELY*COS ) /ROOT  37 
READ  WAVE  PERIOD  FROM  FLEET  NUMERICAL  WEATHER  FACILITY  TAPE  (NOT  USED) 

BUFFER  IN(1»2)  ( ND ( 1 )  ♦  ND( 3969 ) )  38 

lb  IF(UNIT.l)  15.24.18,20  39 

18  WRITE(51.19)  K  40 

19  FORMAT  (41H0  PERIOD  EOF  OR  EOT  ERROR  OCCURRED  ON  K=I3/)  41 
GO  TO  6  42 

20  WRITE(51.21)  K  43 

21  FORMAT  (46H0  PERIOD  PARITY  OR  LENGTH  ERROR  OCCURRED  ON  K=I3/J  44 
M  =  LENGTHF(l)  45 
IF  (M-3969)  22*2^**22  46 

22  WRITE (5 1,23)  M  47 

23  FORMAT  (25H0   PERIOD  BUFFER  LENGTH  =16/)  48 
READ  WAVE  HEIGHT  H  FROM  FLEET  NUMERICAL  WEATHER  FACILITY  TAPE 

24  BUFFER  IN(1.2)  ( ND ( 1 ) ♦ ND( 3969 ) )  49 

25  IF(U.MIT.l)  25»34,28,30  50 

28  WRITE(51»29)  K  51 

29  FORMAT  (41HC  HEIGHT  EOF  OR  EOT  ERROR  OCCURRED  ON  K=I3/)  52 
Nl  =  K  53 
GO  TO  16  54 

30  WRITE(51.31 )  K  55 

31  FORMAT  (46H0  HEIGHT  PARITY  OR  LENGTH  ERROR  OCCURRED  ON  K=I3/)  56 
M  =  LENGTHF(l)  57 
IF  (M-3969)  32*34,32  58 
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32  WRITE(51»33)  M  59 

33  FORMAT  (25H0   HEI6HT  BUFFER  LENGTH  =16/)  60 

34  DO  43  1=1,18  61 
DO  43   J=l,32  62 

IH  =  MD( I+8,J+16)/2048  +  ISCALE  63 

43  XHT(ItJtK)  =  (H  +  0.0)*64.  64 
REWIND  1  65 
WRITE(2»47)  XHT,CSK,SNK  66 

47  FORMAT  (13F9.5)  67 

IF( IOCHECK,2)  44,46  68 

44  WRITE( 51,45  )  69 

45  FORMAT  (45H0   PARITY  ERROR  OCCURRED  ON  XHT,CSK»SNK  WRITE/)         70 

46  END  FILE  2  71 
Nl  =  10  72 

16  WRITE(51,27)  <  (  (XHT(  I ,J,K) ,1=1,17)  ,J  =  1,31,5)  »K  =  1»N1 )  73 

27  FORMAT  (17F7.2/)  74 

6  REWIND  2  75 

STOP  76 

END  77 

END  78 

FINIS  79 

JECUTER.  80 

>78   9781746117461   978   978  2070  4931  4133  3388  2698  2065  1491   978   978ABS  1 
m\    2064  2697  5198  4550  3956  3418  2937  2515  2152  1851  1612  1435  1321  1271ABS  2 
!84  1361  1501  1704  1969  2296  2684  3130  3635  4197  4814  5484  6205  6976  7794ABS  3 
i'64  9841  8966  8144  7378  6669  6021  5437  4917  4465  4082  3769  3528  3359  3262ABS  4 
|!40  3290  3414  3611  3879  4219  4628  5106  5651  6259  6931  766110146  9307  8535ABS  5 
134  7207  6658  5188  5800  5496  5278  5146  5102  5145  5275  5491  5793  6179  6647ABS  6 
.96  7821  8520  929 1 10128 1 1028  1 1987 1 3000  14062 1 7461 1 7461 16477 1552314600 137 14ABS  7 
J68120671131210608  9957  9363  8827  8353  7942  7596  7316  7104  6960  6885  6880ABS  8 
'44  7078  7280  7550  7886  8288  8753  9279  9865 105071 12041 1951 12746 1 567214829ABS  9 
1251326512551118871127  71072310227  9793  9422  9116  8877  8705  8601  8567  8601ABS10 
§5  6877  9116  9422  97931 022 7  1072 31 1277 1 188 7 12551 1 3265 1402614829 1 567 316551 ABS11 
■6 117461 1666315  89  7 1516 7 144751 382  513220 1266412159 117081 1311 1097 31 06941 0475ABS12 
■181 022410 1931022410318 104751069410973 1131 11 1708 12159 12664132201 382 514475ABS13 
67 1^89 7 16663 1746 11 746 11 682 7 1622 11 5645 15 100 14590 141 171 3682 1328 81 2936 12628ABS14 
3 6121 49119791185811785117 60 11 78 51185811 97 912149 123 66 126281 29 361 32881 36 82 ABS15 
171459  15 100 15645 1622 1168271746 11 7461 16989 16539 161 1315 712 15388 1499 114674ABS16 
8 71 41 3 113903 13 71 7 13560 13438 13350 13297 13279 13297 13350 13438 13560 1371 71 3908ABS17 
31 1438 7 1*6 74 149 91 1533815  712 161 13165  39 16  989 17461 17461 17 1501685 5 16577163 17ABS18 
741 585 11564715462 15 2991 5 1561 50 351493 5 1485 7 148 01 14768 14 7 57 1476 8 14801 1485 7ABS19 
351 503515 1561529915462 1564715851 160741631716577 16855171 501746 11 746 11 7312AB520 
72 1704 1169 19 168 06 1670 2 16608 1652 3 16448 16 382 16327 1628 116246 16220 1620 51 6200 ABS21 
0^162  2  162461628 11632 7 16382 16448 16523 16608 16702 16806169 19 17041 171 72 173 12 ABS22 
61174611672215277174611746113665118241746117461  9667  70621746117461  3797ABS23 
W8   9731746117461   978   9781746117461   978   9781746117461   978   978 1 7461ABS24 
561   978   9731746117461   978   9781746117461   978   9781746117461  2520  5878ABS25 
16117461  86421099117461174611304014870174611746116537  ABS26 

5292870528705   629   629  2224   629   629  1612  2637  3700  4797  5927  708525322ORD  1 
I  8 02 76 0  92 8 705 28 70 52 772 126 70 3256 552 45 78 23478 223 57 2 12 172 0063 18898 1772 51 65470RD  2 
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153681*192130221186110713  9580  8467  7377  6312  5275  4271  3302  2369  1478   6290RD 
62^  1408  22-42  3127  4061  5039  6059  7116  8207  9328 104731 1640 12824140201 52240RD 
164321  /bJ9 18  84 1200  3  32 12 112 2  3 7 12  3507246172  5695267392  774328 70528 70 52 777 2267 8 30RD 
2  5  ^2*655  23  3  282  2  3642117  1199  531871817470162161496213714124  78112611006  7  890  3ORD 
7774  b686  5645  4654  3720  2846  2037  1297   629   629   976  1379  1844  2372  29580RD 
3601  4298  5046  5842  6682  7562  8480  9431 10410 1 1415 12440 1348 1 14533 1559 3 166560RD 
1771  bid  7 7  198  1320841 2 18492 2832 2 37882471025597264432 72462800 1287052870 528 1630RD 
2 7 5642b9 11 2620725456246602 382 5229532205021 11820 164 19 190 1820317205 1620 31 5202ORD] 
1420413217122431128810357  9453  8582  7746  6951  6200  5496  4843  4244  3702  32190RD] 
2797  4582  5005  5483  6013  6593  7220  7890  8602  9350 10132 1094311 780 12638135 140RD] 
1440 3 1530  116204 17 106 180041 8893 197692062721464222762 305723806245 1725 18 7258 140RD] 
26 3 942 69242 74 022 782 5 2605 2 256602 522 72 475624248 2 37052 3 130 22 52 52 189 32 123 7205 590RD] 

19  86219149184241768916943162  041545914718139831325812  545118491117010514  98  8  20RD] 
9277  87oZ    8159  7651  7180  6747  6355  8154  8507  8887  9294  97261018 1 1065 71 1 1540RD] 

11 668 12 199  12 74413 30 11 3869 144451 5028 1561 5 1620416792 17379 17962 18538 1910619663ORD] 

20  2^  32  0  7  392  12  5^+2  1750 2222  62  268  123  11 32  3 52023900  242 5322 360220602 17432  141  22  1067ORD] 
20 70 ^2 03^8 1995 7 19 566 19166 18757 18 342 1792 117496 17067 16636 162041 577 11 5340149 110RD: 
1 4486 1^06  5 1365 01 3242 1284 11 2450 12 069 11698 11 340 10995 10664 10 347 10047 1222 2 124420RD; 
1266 7 1289 7  13  133 13372 136161 386414 11 514370 14627 14886 15 147154 10 156741 593816204ORD; 
i6469 167 3 316997 17260 1752 11 778 118038 18292 18543 18791 19035 19275 19510 19740 199660RD; 
20  18 52 726 128 7 0523 70 52 52 2 4236 852 8 70528 70 52 2466 2 1462 2 8 7052 8 70520608 198642 870 50R D; 
28 70 ^2  792 9 19 200 18 5942 55582 33 57 1803 11 75002 1277 19278 16989 16491 173271 53941 59970RD; 
1549913446114591499114464  9395  72161390713309  4876  23181265511926  629  6290RDi 
11C9310H7   629   629  8938  7461   629   629  5521  2814   629  ORD; 

5UBR0UTINE  AP2 

COMMON  YC(10)  ,LR,A,B,CC,H,CK»SK  1 

COMMOf\i/L2/HX»HY»HXX,HXY,HYY  2 
CQMM0N/L4/AX,AY,bX,6Y,CX,CY,AxX»AXY,AYY,BXX,BXY,BYY»CXX,CXY,CYY  3 

Rl    =    SQRTF( ( .041783709356*H-.42321401072)*H+2. 2342337579)  4 

VF    =    -.C28281950577*H    +    17.494735347    -    Rl  5 

R2    =    SQRTF( ( .058458667266*H-.90729449065)*H+6. 4033842487 )  6 

B  =  -0.14520001518*H  +  18.530490911  -  R2  7 

DVFM  =  (-.04T783  709356*H+.21160700536)/R1  8 

DVF  =  DVFM  -  .028281950577  9 

DBM  =  (-.058458667266*H+.45364724533)/R2  10 

DB  =  DBM  -  0.14520001518  11 

D2VF=  (DVFM*DVFM  -  . 041783709356 ) /R 1  12 

02b  =  (DBM  *  DBM  -  . 058458667266 ) /R2  13 

IF  (H-15.)   1,1,2  14 

1  K3  =SQRTF((  . 23341292994*H-3.1096617758)*H+29. 275404601 )  15 
VH  =  -0.25152614353*H  +  21.408836894  -  R3  16 
DVHM  =  (  -.23341292994*H+1.5548308879)/R3  17 
DVH  =  DVHM  -  0.25152614353  18 
D2VH=  (DVHM*DVHM  -  . 2 3341292994 )/ R3  19 
GO  TO  3  20 

2  R4  =SQRTF((  . 14668786198*H-6. 8828319323 ) *H+1 05 . 12448592 )  21 
VH  =  -0.36970234218*H  +  11.346369501  +  R4  22 
DVHM=  (  0.14668786198*H  -  3. 4414 1596615 ) /R4  23 
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OVH  =  DVHM  -  0.36970234218 

D2VH=  (-DVHM*DVHM  +  . 146687861 98 ) /R4 

A  =  (VF+VH)*.5 
CC  =  A-VH 
DA  =(DVF+DVH)*.5 
DC  =DA-DVH 
AX  =  DA*HX 
AY  =  DA*HY 
bX  =  DB*HX 
BY  =  DB*HY 
CX  =  DC*HX 
CY  =  DC*HY 
IF  (LR)   4*6»4 
XX  =  HX*HX 
XY  =  HX*HY 
YY  =  HY*HY 
D2A  =( D2VF+D2VH)*.5 
D2C  =  D2A  -  D2VH 
AXX=  DA*HXX  +  D2A*XX 

DA*HXY  +  D2A*XY 

DA*HYY  +  D2A*YY 

DS*HXX  +  D2B*XX 

DB*HXY  +  D2B*XY 

D8*HYY  +  D26*YY 

DC*NXX  +  D2C*XX 

DC*HXY  +  D2C*XY 

DC*HYY  +  D2C*YY 


AXY  = 
AYY  = 
bXX  = 
8XY  = 
OYY  = 

cxx= 

CXY  = 
CYY  = 
RETURN 
END 


24 
25 
26 
27 
28 
29 
30 
31 
32 
33 
34 
35 
36 
37 
38 
39 
40 
41 
42 
43 
44 
45 
46 
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48 
49 
50 
51 
52 


List  of  cards  requiring  changes  if  the  geometric  dimensions  and/or 
origin,  or  time  dimension  of  the  XHT ,  CSK  and  SNK  arrays  are  changed: 

VC2AP3:  2,3,4,7,15,27,44,94  to  97,113,114,117,118,146,147,162,163, 
168,169,212,213,291,292,306,307,334,335,339  to  341 

TERP:    5,6,38,40,50,55,121,133 

ANGLE:   1  to  4 

VTJERIV:  1  ,8,9 

AP3  (and  substitutableAP2) :  1 

TAPE:    2,3,8,15,27  to  30,32,61  to  63,72,73 

Some  unnecessary  remnants  of  experimentation  remain  to  be  cleaned 
up  as  follows : 

TERP:         When  LR=0  only  H,  CK  and  SK  need  be  computed  before 
returning  to  the  main  VC2AP3  program 

AP3  and  AP2:  Vhen  LR=0  only  A,  B  and  CC  need  be  computed  before 
returning  to  the  main  VC2AP3  program 


VDERIV: 


When  LR=0  only  CAPV  need  be  computed  before  returning 
to  the  main  VC2AP3  program 
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G 1 ossary  (For  o 

X(900)  ,Y(900)  = 

1T(12) 
C(4) 

AK(4,10) 

DY(IO) 

D(20) 

LR 

A,B,CC 
H  ,  CK  ,  SK 

XHT , CSK , SNK 

COST,SINT 

V 

CAPV 

CAPVX  to  CAPVYY: 

VPBVX  to  VPBVY= 

DIV 

RBV 

YVARS(10) 


missions  consult  pages  9  and  10) 
VC2AP3 

Storage  for  map  grid  data,  and  temporary  storage 

for  output  data 

Title  for  map  produced  by  DRAW  subroutine 

Runge-Kutta  integration  weighting  factors 

Data  used  in  finding  four  YC(10)  intermediate 

ordinates  in  one  Runge-Kutta  integration  step 

Time  derivatives  in  the  10  differential  equations 

Store  used  in  integrating  eqs.  for  6x,6y,6?,6£ 

Index  to  denote  a  geodesic  track  (LR=0)  or  a 

minimal-time  track  (LR=1 ) 

Parameters  for  elliptical  polar  velocity  figure 

Wave  height  and  wave  direction  cosine  cosK  and 

sine  sinK  relative  to  Oxy  axes 

Wave  field  arrays  with  18x32x10  or  5760  elements 

cos0,  sinS  of  the  polar  velocity  diagram 

Speed  in  knots  on  earth's  surface 

Speed  in  grid  units  per  day  in  stereographic  plane 


XK,XLG,XLT 


RSTEP 

C0SLG1 ,SINLG1 

C0SLG2,SINLG2 

C0SLT1 ,SINLT1 

C0SLT2,SINLT2 

EL , EM , EN 

X1  ,Y1  ,X2,Y2 

S1  2 

ARC  (Card  81  ) 

ALF 

C0SA,SINA 

XFIN,YFIN 

STEP 

TAU 

S 

TVAR 

TC 

N1 

N2 

N3 


^Partial  derivatives  of 
Partial  derivatives  of 

(V  /V)2-(V  /V)  +1 

P      O    Pi    P 

[(vp/v)M* 


CAPV 
V  /V 


X1  ,M.  ,\~  * M-p  ,x,y,  geodesic  arc  length  or  6x,6y,6^, 


1 


and  6£  at  end  of  each  integration  step 

Vave  direction  from  north,  longitude  and  latitude 

produced  by  ANGLE  subroutine  for  printed  output 

tabulation  by  statements  72  to  74  of  VC2AP3 

Floating-point  value  of  NSTEP 

Cosine  and  sine  of  longitude  of  initial  point 

Cosine  and  sine  of  longitude  of  terminal  point 

Cosine  and  sine  of  latitude  of  initial  point 

Cosine  and  sine  of  latitude  of  terminal  point 

Direction  cosines  (later  normalized)  of  normal 

vector  to  great  circle  track  plane 

Coordinates  of  initial  and  terminal  points  of 

track  relative  to  North  Pole  in  grid  units 

Straight-line  distance  from  X1 ,Y1  to  X2,Y2 

Great-circle   distance  from  X1 ,Y1  to  X2,Y2 

a  =  Departure  angle  of  ship  track  at  initial 

point  relative  to  Oxy  axes  as  in  Figure  l(approx.) 

Cosine  and  sine  of  ALF 

Coordinates  of  desired  terminal  point  relative  Oxy 

Time  integration  step  in  units  of  days 

Time  in  days  from  beginning  of  track 

Arc  length  on  geodesic  track 

Time  in  days  from  first  member  of  time  series 

Intermediate  values  of  TVAR  in  one  integration  step 

Index  for  number  of  integration  steps 

Index  used  in  plotting  a  geodesic  track 

Index  used  to  indicate  whether  a  complete 

geodesic  track  can  be  computed 
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XN1 
FSTEP 

XLAM 
XMU 

CLAM 

COSP,SINP 
FAC1 
DET 

EMFI 

XDOT,YDOT 
DIFT,DIFA 
NK 


Time  in  days  to  complete  geodesic  track  or 
any  time-extremal  iterated  track 
Floating-point  value  of  N1 

Final  Runge-Kutta  process  time  increment  on  a 
time-extremal  track,  with  STEP1  used  otherwise 
>  =  X-,  cosa  + 
|i  =  u1  cosa  +  | 
2.2" 


XpSina 
UpSina 


=  Direction  cosines  of  CAPV  relative  to  Oxy  axes 
to  FAC4   =  Temporary  storage 

=  \^"?    ~      ^^  1 

=  Ratio  of  a  differential  distance  in  stereographic 

plane  to  corresponding  distance  on  earth's  surface 
=  x  and  y 
=  AT  and  6a 
=  Index  for  daily  position  point  plot 

TERP 


C(4) 


L 

TT 

M,N 

XX,  YY 

P(1 )  to  P(4) 

Q(D  to  Q(4) 

PX,QY,PXX,QYY 

H 

CK,SK 

HX  to  SKYY 

DKX  to  DKXY 


VF 

VH 

A,B,CC 

AX  to  CYY 

SIN2T,C0S2T 

EMFIX,EMFIY 

VBR 


VPBVP 


Weighting  factors  for  interpolation  in  the  time 

dimension  between  K4  ordinates 

Index  to  pick  out  member  of  wave  field  time  series 

t  of  Eq.(29)  of  NPS  Tech. Report/Res .Paper  No.  46 

Indicies  to  pick  out  x,y  grid  point  data 

x,y  of  equation  (18)  of  this  report 

P,  to  P.  of  equation  (19)  of  this  report 

Elements  of  P'(y)  matrix  of  Eq.(l8)  of  this  report 

Partial  derivatives  of  P(1  to  4)  and  Q ( 1  to  4) 

Interpolated  wave  height  in  feet 

Interpolated  (and  later  normalized)  cosK  and  sinK 

Partial  derivatives  of  H,  cosK  and  sinK 

Partial  derivatives  of  the  wave  direction  angle  K 

AP3 ,  AP2  and  VDERIV 

Ship  speed  in  knots  in  presence  of  following  waves 

Ship  speed  in  knots  in  presence  of  head  waves 

Parameters  for  elliptical  polar  velocity  diagram 

Partial  derivatives  of  A,  B  and  CC 

sin  20 ,  cos  29 

Partial  derivatives  of  EMFI 

Ship  speed  V  in  knots,  multiplied  by  24  hours,  and 

divided  by  stereographic  plane  mesh  size  of  205.585 

nautical  miles  at  60N  latitude 

(V  /V) 
P    P 
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